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Summary 

A summary of the work performed under NASA grant NCC3-605 is presented. More details can 
be found in the cited references. This grant led to the development of relatively faster aeroelastic 
analyses methods for predicting flutter and forced response in fans, compressors, and turbines 
using computational fluid dynamic (CFD) methods. 

Linearized Euler Solvers 

Previous aeroelastic analyses have used the unsteady aerodynamic equations based on uniform 
flow. These methods neglect the effect of steady aerodynamic loading on unsteady aerodynamic 
forces. Steady loading is due to angle of attack and airfoil shape effects. These methods require 
the least amount of computational time to perform an aeroelastic analysis. Advanced methods 
based on computational fluid dynamic (CFD) methods have been developed to include steady 
loading effects in the aeroelastic analysis. However, these methods, usually based on time 
marching, require an enormous amount of computational time for the aeroelastic analysis. 
Methods based on the linearized, unsteady, aerodynamic equations combine the benefits of the 
above two methods. They are faster and include steady loading effects. This is achieved by 
linearizing the non-linear unsteady equations using a non-linear steady aerodynamic solution to 
obtain linear, unsteady, aerodynamic equations. These equations can be solved in the frequency 
domain, by assuming harmonic motion for the unsteady oscillations. The resulting analysis code 
is faster and includes steady loading effects. 

In this grant period, aeroelastic analysis codes have been developed using a two-dimensional 
(2D) linearized Euler solver, LINFLUX-2D, and a three-dimensional (3D) linearized Euler 
solver LINFLUX-3D. The difference between them is in using different structural models. 

Also, during this grant period, a version of the LINFLUX-2D code and a version of the 
LINFLUX-3D code were implemented on an SGI machine, and the unsteady aerodynamic 
calculations were validated with published results. The execution of LINFLUX-3D required 
learning and running a steady aerodynamic code, TURBO-AE on SGI machines. This in turn 
required the development of an interface code to link the steady solution from the TURBO-AE 
code to the LINFLUX-3D code. 

Two-Dimensional Aeroelastic Analysis code MISER- LE 

A typical section structural model was used to develop the aeroelastic code, MISER-LE. The 
unsteady aerodynamic forces are obtained from the based on linearized Euler solver, LINFLUX- 
2D. This structural model has two degrees of freedom, bending and torsion. The governing 
aeroelastic equations are solved in the frequency domain for each interblade phase angle. 
Cascades with 9 and 12 blades were considered in the analysis. The unsteady aerodynamic 
forces were obtained from linear theory, where steady loading effects were neglected, and from 
LINFLUX-2D where these effects were included. The results are presented in Ref. 1. An 
extended version of this paper, showing an extension to forced response and correlation with 
previously published data is given Ref. 2. 
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Quasi-3D Aeroelastic Analysis code, ASTROP2-LE 

In this analysis, a three-dimensional structural model is used for the aeroelastic analysis. The 
unsteady aerodynamic forces are obtained from LINFLUX-2D. The structural and the unsteady 
aerodynamic models are combined using strip theory. An existing code, ASTROP2 was taken 
and updated to include forced response, and the unsteady aerodynamic solution from LINFLUX- 
2D. The governing equations are solved in the frequency domain. The details of the resulting 
code ASTROP2-LE and the results for a tuned cascade aeroelastic analysis were presented in 
Ref. 3. The study was later extended to include mistuning effects and is being reviewed for 
NASA TM publication, Ref. 4. 

Three-Dimensional Aeroelastic Analysis code, LINFLUX-AE 

In this effort, an aeroelastic system, LINFLUX-AE, is being developed for the aeroelastic 
analysis of three-dimensional structures. A three-dimensional structural model was combined 
with the three-dimensional, unsteady, aerodynamic model LINFLUX-3D. A normal mode 
approach combined with the frequency domain solution method is used for aeroelastic analysis. 
Modules were developed to interpolate structural mode shapes on aerodynamic grids, calculation 
of generalized forces and flutter eigenvalues. The modules were checked with known results. 
The details of the preliminary version of the LINFLUX-AE code were presented in Ref. 5. The 
paper presents flutter calculations for a helical fan with flat plate geometry, and for a real fan, E- 
cubed fan. Flutter eigenvalues and work done per cycle were compared with those obtained 
using TURBO- AE. Further validation and extension to forced response are under progress. 
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ABSTRACT 

Flutter analysis is presented for a cascade of 
blades in subsonic and transonic flow. The 
structural model for each blade is a typical section 
with bending and torsion degrees of freedom. 
The unsteady aerodynamic forces are obtained by 
solving unsteady linearized Euler equations. The 
unsteady linearized equations are obtained by 
linearizing the unsteady non-linear equations 
about the steady flow. The predicted unsteady 
aerodynamic forces include the effect of steady 
aerodynamic loading due to airfoil shape, 
thickness and angle of attack. The unsteady 
aerodynamic and aeroelastic equations are solved 
in the frequency domain. The present aeroelastic 
solver showed good correlation with published 
results. Further improvements are required to use 
the unsteady aerodynamic model in a design 
cycle. 

INTRODUCTION 

The aeroelastic research program at NASA Lewis 
Research Center is focused on flutter (unstalled, 
stalled, and whirl), and forced response analysis of 
propulsion components. An overview of this 
research was presented in Ref. 1 . The review 
showed that a range of aerodynamic and structural 
models have been used to obtain the aeroelastic 
equations. Both, time and frequency domain 
methods have been used to obtain unsteady 
aerodynamic forces and to solve the aeroelastic 
equations. It was noted that time domain methods 
require large computational time compared to 
frequency domain methods, and should only be 
used when non-linearities are expected, and may 
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be used in the final design. 

Two approaches have been used in obtaining the 
unsteady aerodynamic forces using frequency 
domain methods. In the first approach, Refs. 2-6, 
the unsteady aerodynamic equations are 
linearized about a uniform flow, there by 
neglecting the effects of airfoil shape, thickness, 
and angle of attack. The unsteady aerodynamic 
models developed in Refs. 3-6 have been used to 
study the flutter and forced response analysis of a 
compressor cascade using a typical section 
structural model, in Refs. 7-8. However, methods 
developed by this approach are restricted to 
shock-free flows through lightly-loaded blade 
rows. In the second approach, Ref. 9, the 
unsteady flow is regarded as a small amplitude 
perturbation about a non-uniform steady flow. The 
unsteady non-linear aerodynamic equations are 
linearized about the non-uniform steady flow, 
resulting in variable coefficient linear unsteady 
aerodynamic equations, which include the effects 
of steady aerodynamic loading due to airfoil shape, 
thickness and angle of attack. 

Following the second approach, Refs. 10-11 
developed a nonlinear steady and linear unsteady 
aerodynamic model based on the potential 
equation. This unsteady aerodynamic model was 
used to study the effect of steady aerodynamic 
loading on flutter stability using a typical section 
structural model in Ref. 12. However, the 
formulation based on the potential equation 
requires corrections for entropy and flow rotation. 
The Euler equations can be used to correctly 
model rotational and entropy effects associated 
with strong shocks. Unsteady linearized Euler 
aerodynamic models that include the effect of 
steady aerodynamic loading were developed in 
Refs. 13-15. 

Recently, a two dimensional linearized Euler code 
named LINFLX2D was developed in Ref. 16 under 
a NASA contract. This code is based on the non- 
linear Euler solver, NPHASE, developed in Ref. 
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17. The objectives of the present study are (1) to 
couple the LINFLX2D code with the aeroelastic 
code MISER (Refs. 7,8 ), (2) validate the 
aeroelastic predictions, and (3) evaluate LINFLX2D 
code for its usefulness in a aeroelastic design 
cycle. The MISER code at present calculates the 
aeroelastic stability using the unsteady 
aerodynamic models based on classical linear 
theories (theories based on flat plate geometry 
formulation). MISER code is selected because the 
typical section structural model with pitching and 
plunging degrees of freedom is the basis for 
aeroelastic formulation. This is in line with the 
unsteady aerodynamic model of LINFLX2D. In the 
present paper flutter analysis of a cascade of 
blades is presented using the frequency domain 
approach for a two dimensional cascade model. 
Two degrees of freedom, plunging and pitching, 
are considered in the analysis. Brief descriptions 
of the formulation and method of analysis are given 
in the next section, followed by results and 
discussion. 


FORMULATION 

The aerodynamic model and the aeroelastic 
formulation are described in this section. 

Aerodynamic model 

Non-linear Steady Euler Solver. NPHASE 

The steady aerodynamic model is based on the 
unsteady, two dimensional, Euler equations. The 
equations in conservative differential form are 
solved in a time dependent body-fitted curvilinear 
reference frame. This transformation process and 
the ensuing numerical method are presented in 
detail in Ref. 17. The equations are discretized 
and solved using a finite volume method using a 
TVD scheme. The steady solutions presented 
herein are obtained using the basic implicit 
scheme developed in Ref. 17, which is third order 
accurate spatially and second order accurate in 
time. 

Linear Unsteady Euler Solver. LINFLX2D 

To obtain the linearized unsteady Euler equations, 
Ref. 16, the independent variables in the 
unsteady non-linear Euler equations are 
expanded in an asymptotic series of the form 


U = U(x) + u(x(x,t),t) + higher order terms (1 ) 

Assuming the unsteady excitations are harmonic in 
time, and with the first order variable to be 
represented as complex valued, the above 
equation can be written as 

U=U(x) + Re[u(x)exp(iwt)] + higher order terms (2) 

Here, the term U(x) is of order one and the second 
term is of the order e. 

Substituting the expansion of Eqs. 1 and 2 in the 
nonlinear unsteady Euler equations, and equating 
terms of like power in e, and neglecting terms of 
second and higher order in e, nonlinear steady 
equations and linear variable coefficient unsteady 
equations are obtained. 

For harmonic blade motions with constant phase 
angle between adjacent blades (interblade phase 
angle), the values of interblade phase angle (s) 
that can occur are given as (Ref. 18). 

o,-2nr!N r = 0 , 1 , 2 , ... ,N - 1 (3) 

N is the number of blades in the cascade. In a time 
domain approach, the number of blocks required 
depends on the inter blade phase angle, and small 
phase angles may require large number of blocks 
to calculate the unsteady aerodynamic forces. 
However, with the linear approach, the periodic 
conditions are applied on a single extended blade 
passage region i.e. a region of angular pitch, 

0 = 2* n/N (4) 

In solving the linear unsteady equations, the 
independent variables are regarded as pseudo 
time dependent. This allows solutions to be 
determined using conventional time -marching 
algorithms to converge the steady and the 
complex amplitudes of the unsteady conservation 
variables to their steady state values. For more 
details, see reference 16. 

Grid and Boundary Conditions 

The flow equations are solved on a passage- 
centered H-grid. Two grid blocks are shown in Fig. 
2 for calrity. Within a typical grid block, the lower 
computational boundary contains the upper 
surface of one blade in the cascade, while the 
upper computational boundary contains the lower 
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surface of the adjacent blade. Periodic boundaries 
in the blade-to-blade direction extend upstream 
and downstream from the blade surfaces. The 
inlet boundary corresponds to the left 
computational boundary and the outflow 
corresponds to the right boundary. The airfoil 
upper and lower surfaces are located along lines B- 
C and F-G, respectively where solid wall boundary 
conditions are employed. Lines A-E and D-H 
represent inflow and outflow boundaries, 
respectively. For subsonic flow these conditions 
are set using characteristic variable boundary 
conditions. The procedure used is to fix the 
incoming flow incidence angle and adjust the back 
pressure (uniform across D-H) until the average 
Mach number along the inflow boundary (A-E) 
matches some specified value. For supersonic 
flow, the inlet conditions are assumed to be 
uniform by specifying the flow density, velocity, 
flow angle, and pressure. The outflow variables 
are extrapolated from the interior by using simple 
first-order model. Periodicity is imposed between 
lines A-B and E-F, and lines C-D and G-H. 


Aeroelastic model 

For completeness, the aeroelastic formulation in 
the frequency domain ( Refs. 7,8) is presented in 
this section. The analysis is presented for a 
generally mistuned cascade in which each blade 
may have different structural properties. The 
analysis for the special case of a tuned cascade in 
which all blades are identical is presented in the 
subsequent section. The approach followed 
assumes that the structure is vibrating in an 
aeroelastic mode (interblade phase angle mode) 
with a motion that is a harmonic function of time. 
The frequency of oscillation is permitted to take on 
complex values thus allowing decaying-, growing- 
or constant-amplitude oscillations. The 
aerodynamic forces corresponding to constant- 
amplitude harmonic oscillations are inserted Into" 
the equations of motion to formulate a complex 
eigenvalue problem. The eigenvalues are 
generally complex quantities, and therefore a 
complex frequency is obtained. The real part of the 
complex frequency represents the damping ratio 
and thus its sign determines whether the motion is 
decaying or growing; the imaginary part represents 
the damped frequency of oscillation. 


Aeroelastic analysis for mistuned cascades 

The equations of motion for the typical section 
(see Fig. 1) with structural damping can be written 
in matrix form for the s" 7 blade as: 

[Ms 1 ,9s , + [C s ] fls } + \K s }q s \- ( s = {/as} + \f B s (5) 

where {f a ^ ) denotes motion-dependent 

aerodynamic loads and { f a ) denotes motion- 
independent aerodynamic loads. 

For the two degrees of freedom considered here, 
equation 5 can be rewritten as 

1 *as j h s lb | ICQjis^hs 0 jh s lb\ 

x as r as_ t a $ ( 0 2r^ffl os £ cts | tts j 

+ (Ohs 0 fh s /b j _ j fhg/mgb \ 

. 0 rlsO^as l a s j \fas/m s b 2 1(6) 

where h is the plunging (bending) displacement 
normal to the chord, a is the pitching (torsion) 
displacement, x a the distance between the elastic 
axis and center of mass in semi-chord units; r a is 

the radius of gyration about the elastic axis in semi- 
chord units; £/j and are the damping ratios; b is 
the airfoil semi-chord; is the uncoupled natural 
frequency for bending; co a is the uncoupled 
natural frequency for torsion; ff j and f a are the 
aerodynamic loads and s varies between 0 and N- 
1 . 

It is assumed that the motion of the blades is 
harmonic in time with a frequency co and is given by 

'v * 1 

I I I I „0 l “" r I 

Note that the motion has been represented as the 
sum of contributions from each interblade phase 
angle mode in which each blade has an amplitude 
h af /b, a ar and the phase angle between adjacent 
blades is given by Eq. 3. 

The corresponding aerodynamic forces can be 
written in terms of the complex-valued unsteady 
aerodynamic coefficients Ihh , I ah > !ha< laa . >wh< 
and ! wa . 
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fhs/l^l J} \ 

i/jf + lj,a- Onr 


U r K r ^ + Lar a o r 


where {Y} consists of the displacement amplitudes 
2 N-l corresponding to the rth interblade phase angle 

£ £ v | e '«ty& s and [E] is the transfer matrix, the elements of which 
Ps Imgr are given by 


( 8 ) 

2 

where /i s = m s ' ftp°ob is the mass ratio of the 
blade. 

Using Eq. 7 and Eq. 8, Eq. 6 can be written as 


E(s,r)=e 2nsri ,N 

where s is the blade counter and r is the interblade 
phase angle counter. 

Using this relation, we obtain: 


-[M s }\ h ° slb \e^ + A [iT s ]( Aos/ 6 L^ 
\ / \ OC OS j 


N-l 

r=0 

where 


\hjb\ ■ 


«or 


e l « s e m + 


N-l 

X |AD r je i »*e te 

r= 0 (9) 


[M s ) = p i 
[K s ] = 


1 

r\ 


{(QhsICQo) (1+2 i£hs) 
0 


0 

r as(®ast®o) (l+2j£ os ) 


[Ar] = 


ihhr i-har 
lahr l aar i 

[AD r | = I j Whr I 
I Lwar I 


A =(co 0 loo) 


<do = reference frequency and the damping terms 
are approximated as, 


-Iffmirm [Ej'pp](r|= [Ai|y|+fAD; (12) 

where 

[M] is the mass matrix for all blades with diagonal 
elements consisting of mass matrix of each blade. 
Similarly [K] is the stiffness matrix and [A] is the 
aerodynamic matrix for all blades, {AD} is the 
aerodynamic forces due to wake on all the blades. 

Finally, after rearranging, the equations can be 
written as: 

[[P]-A[Q]]jF) = -(AD} (13) 

where 

[p] = [[£r’[AfFMA]' 

[<J ] = [£]-■ [xp] 

For a stability calculation (flutter), the motion- 
independent forces {AD} are set to zero and the 
eigenvalue problem is obtained in the standard 
form: 


liCQCOhsChs ~ AicOfig C,h s 

Air as COCOas — 2 iVas (Oas^as (-jo) 

To proceed further, the equations for all the N 
blades on the disk must be considered. For the 
assumed harmonic blade motion, the 
displacements {X} of all blades can be written as a 
sum of contributions from all interblade phase 
angles as 

[X]e iat = [£] {Y ) e im (11) 


[[p]-a[q]](y) = (o} (14) 

The solution of the above eigenvalue problem (14) 
results in 2N complex eigenvalues of the form 

i SL = —L — = ji + iv 
a ° IT (15) 

The real part of the eigenvalue ( pi ) represents 
the damping ratio, and the imaginary part ( y ) 
represents the damped frequency; flutter occurs if 

fl>0 for any of the eigenvalues. 
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The aeroelastic response of the blades induced by 
wakes is calculated from equation (13) as 


where 



[PMMf’jAD) 


(16) 


The amplitude of each blade is obtained by 
substituting equation (16) into equation (11). 


[Pr) = 


/f + ^hhr 


HWQa) (1+2 »&) 
lirl (l+2if 0 ) 


(1% c + I'hcn' 

2 

fj ((Qh/cOtx) (1+2 ift) 

+ (g 07' 

M ^'a (l+2iCa) 


Aeroelastic analysis for a tuned cascade 

For a tuned cascade (or rotor), in which all the 
blades are identical, the foregoing analysis can be 
simplified considerably. In this case, the 
aeroelastic modes consist of individual blades 
vibrating with equal amplitudes with a fixed 
interblade phase angle between adjacent blades. 
Hence, for this problem, the motion of the typical 
blade is written as 



Since the blades are identical, the same equation 
is obtained for each blade. Thus, no additional 
information can be obtained by assembling the 
equations for all the blades on the disk as was 
done for the general mistuned system. Instead, 
equation (19) is solved for A/different values of the 
interblade phase angle given by equation (3). As 
before, the equations for the forced-response 
problem are obtained by setting the motion- 
dependent forces to zero; the equations for the 
flutter problem are obtained by setting the motion- 
independent forces to zero. 


For the stability calculation, the equation can be 
simplified as 


PVM[/j]{y} = [o 


(20) 


where the subscript 's' identifying the blade has 
been dropped and the reference frequency u> 0 

has been chosen to be equal to the torsional 
frequency co a . 

The solution of the above eigenvalue problem 
results in two complex eigenvalues of the form 

fd + i v , and flutter occurs if fJ. sS 0 . For the tuned 
cascade, the stability of each phase angle mode is 
examined separately. Hence, the interblade 
phase angle is fixed at one of the values given by 
equation (3), and the 2X2 eigenvalue problem is 
solved. The value of interblade phase angle is 
then changed, and the procedure is repeated for 
each of the N permissible values. The critical 
phase angle is identified as the one which results 
in the lowest flutter speed. 

Stability calculation 

The aerodynamic coefficients are calculated 
before the eigenvalue problem can be set up and 
solved. Since the unsteady aerodynamic 
coefficients depend on the frequency of 
oscillation, it is necessary to assume a frequency co 
(reduced frequency of blade vibration based on 
chord, k c ) in advance to be able to calculate the 

aerodynamic coefficients. In actual calculations, 
the aerodynamic coefficients are functions of inlet 
Mach number Moo , and interblade phase angle 
Of , in addition to cascade geometric parameters. 
In the present study, for a given inlet Mach 
number, and the reduced frequency is varied until 
the real part of one of the eigenvalues M 
becomes zero while the real parts of the remaining 
eigenvalues are negative. The assumed flutter- 
reduced frequency k C f and the calculated flutter 

frequency Vf are both based on (of . Thus, 
these two can be combined to eliminate cof and the 
flutter speed is obtained, namely, 
Vf = VfC C0 o / k c f . Since the inlet Mach 
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number is known, this flutter speed gives the inlet 
condition (speed of sound, a «, ) at which the 

cascade will be neutrally stable for the given Mach 
number. This procedure can be repeated to 
obtain a plot of flutter speed versus Mach number. 
Knowing the operating conditions, it is possible to 
determine whether flutter will occur within the 
operating region and if so, the Mach number and 
frequency at flutter. It should be noted that the 
analysis of a mistuned cascade requires the 
solution of the equations for all phase angles at 
one time for a given reduced frequency. 


RESULTS AND DISCUSSION 


The analysis methods presented earlier are now 
applied to investigate the behavior of cascades 
oscillating in pitch and plunge. First, unsteady 
pressure distributions and loading for selected 
cascade configurations are presented and 

compared with published results to help validate 
the present linearized unsteady Euler solver. 
Finally, flutter predictions with and without 
mistuning are presented. For the calculation of 
unsteady aerodynamic coefficients, first NPHASE, 
the non-linear Euler solver is used to generate the 
steady solution. This steady solution is then used 
as an input to the LINFLX2D to calculate the 
unsteady aerodynamic forces due to small 
unsteady perturbations in pitch and plunge. 

Code Validation 

Calculations were made for a flat plate cascade, 
and a cascade comprising of airfoils designated as 
the tenth standard configuration (CIO), Ref. 16. 
Both the cascades have a stagger angle of 45 
degrees, and a gap to chord ratio (s/c) of unity. 
The tenth standard configuration airfoils are 
constructed by superposing the thickness 
distribution of a modified NACA 5506 airfoil on a 
circular-arc camber line. See Ref. 16 for more 
details. Two Mach numbers are considered, 
M=0.7 and M=0.8 for a reduced frequency of 
oscillation based on semichord (kj of 0.5. For 
M=0.8, the flow is transonic with a shock at 25% 
chord from the leading edge of the airfoil. The 
steady angle of attack is zero for the flat plate 
geometry, and ten degrees for the CIO airfoil for 
subsonic flow (M=0.7), and thirteen degrees for 
transonic flow (M=0. 8). A H-grid of 141x41 grid is 
used for the study. There are 80 points on the 
airfoil, and the inlet boundary is 5 chords from the 


airfoil leading edge and the exit boundary is 1 0 
chords from the airfoil trailing edge. Some 
calculations were also done with a 151x41 grid for 
which the inlet and exit boundaries are located at 
one chord from the leading and trailing edges 
respectively. There are 55 points on the airfoil for 
this grid. 

Unsteady Aerodynamic Pressures 

Subsonic Inflow 

Figure 3a shows the unsteady pressure difference 
distribution obtained for the fiat plate geometry 
operating at M=0.7. The interblade phase angle, 
a, is 180 degrees, the reduced frequency of 
oscillation based on semichord is 0.5, and the 
Mach number of 0.7. The blades are oscillating in 
pitch about the midchord. For the inlet Mach 
number of 0.7, the flow is shock free. Figure 3a 
shows predictions form linear theory (Ref. 3) and 
from the present linearized Euler code, LINFLX2D. 
The predictions from LINFLX2D correlate well with 
the linear theory results. 

Figure 3b shows the unsteady pressure 
distribution obtained for the standard cascade 
configuration, Cl 0. The flow to the cascade is at 
10 degrees angle of attack. Again, the interblade 
phase angle is 180 degrees, the reduced 
frequency of oscillation based on semichord is 0.5, 
and the Mach number of 0.7. The blades are 
oscillating in pitch about the midchord. For the inlet 
Mach number of 0.7, the flow is shock free. Figure 
3b shows predictions form linear theory (Ref. 3), 
nonlinear Euler (Ref. 17) and from the present 
linearized Euler code, LINFLX2D. The predictions 
from LINFLX2D correlate well with the nonlinear 
Euler results. As expected, they show 
discrepancy with linear theory since the effects of 
airfoil geometry, and angle of attack are not 
included in the linear theory. 

Transonic Inflow 

For the standard configuration, CIO, and for an 
inlet Mach number of 0.8, the flow is transonic with 
a normal shock occurring in each blade passage. 
The steady angle of attack is 13 degrees. The 
steady Mach number distribution for this flow 
obtained from NPHASE is shown in Fig. 4a which 
shows a normal shock at about 28% of the chord. 
Results from the full potential solver, Ref. 19, are 
also included for comparison. Both results agree 
well, except that the nonlinear steady Euler solver, 
NPHASE, predicts the shock location slightly 
downstream of that predicted by the full potential 
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solver of Ref. 19. A shock fitting procedure was 
used in Ref. 19, whereas the shock is captured 
naturally in the present solver. The unsteady 
pressure distribution is shown in Fig. 4b, along 
with a comparison with linear theory (Ref. 3) and 
from the unsteady nonlinear Euler solver. The 
unsteady results are for pitching about mid-chord 
with kb = 0.5, a =180° and amplitude of 
oscillation, a 0 ,= 2°. As expected, linear theory 
does not show any shock, indicating that unsteady 
analysis based on linear theory will not be accurate 
for cascade flutter analysis in transonic flow. 

Similar results were obtained for plunging motion 
for both subsonic and transonic inflows. 

Unsteady Aerodynamic Force coefficients 

It is the unsteady aerodynamic force coefficients, 
unsteady lift and moment that are obtained by 
integrating unsteady pressure differences, that are 
used in the flutter and forced response prediction, 
Table 1 shows the unsteady aerodynamic force 
coefficients obtained using the Smith code, Ref. 
3, and the present LINFLX2D code. Calculations 
were also done with the non-linear Euler solver of 
Ref. 19. In table 1, the first column shows the 
geometry of the airfoil studied, the fourth and fifth 
columns show the force components for pitching 
motion and plunging motion respectively. The 
symbols in the second column are cl for lift 
coefficient, and cm for moment coefficient, with 
moment taken about the mid chord (x0=0.5). The 
first row shows the calculations for the flat plate 
geometry in subsonic flow, the second row for the 
airfoil in subsonic flow, and the third row for the 
airfoil in transonic flow. The notation indicating the 
analysis method is given at the bottom of the table. 
The geometry and unsteady conditions are given 
on the top of the table. 

A comparison of the force coefficients for flat plate 
geometry in the first row shows that the prediction 
by LINFLX2D are very close to those given by 
Smith code. For a flat plate geometry in subsonic 
flow, the Smith code gives the exact answer. It is 
to be noted that the flat plate is unloaded and 
shock free. The values of the force coefficients in 
the second row show some difference from the 
Smith code, since in LINFLX2D the effect of airfoil 
and angle of attack are included. The force 
coefficients in the third row include the effect of 
shock and its motion in addition to the effects of 
airfoil shape and angle of attack. The force 
coefficients obtained with LINFLX2D code 
compare well with those obtained from a non-linear 
unsteady Euler solution, indicating that LINFLX2D 


code can be used to analyze cascades in transonic 
flow. It is to be noted that with 151x41 grid, 
solutions for transonic flow could not be obtained 
with LINFLX2D code. It was suggested that for the 
LINFLX2D code to work properly, it is better if the 
inlet and exit boundaries are at least 5 chords from 
the leading and trailing edges respectively (Ref. 
20 ). 

Flutter Calculations 

As mentioned before, given the unsteady 
aerodynamic force coefficients, the MISER code 
calculates the aeroelastic stability of tuned and 
mistuned cascades oscillating in pitch and plunge. 
The aeroelastic stability results from MISER are 
presented as a root locus plot, with real part of the 
eigenvalue indicating damping and the imaginary 
value indicating frequency for all possible 
interblade phase angles. In the present study 
LINFLX2D is loosely coupled with MISER (Ref. 7) 
to calculate flutter of mistuned cascades. First a 
steady solution is obtained from NPHASE for the 
required cascade geometry and flow conditions. 
With the steady solution as input, LINFLX2D is run 
for given a reduced frequency, interblade phase 
angle, and mode of interest (pitching or plunging). 
The unsteady aerodynamic coefficients obtained 
from LINFLX2D are then stored in a data base, and 
later used in MISER to calculate the stability 
characteristics. 

The structural properties are the mass ratio, m, is 
258.50, ratio of bending to torsion frequencies w/j 
iw a is 0.357, the radius of gyration about the 
elastic axis in semi-chord units, r a ,is 0.5774, with 
no structural damping. The elastic axis is at 
midchord i.e. ah = 0.0. The distance between the 
elastic axis and center gravity, x a> is zero, 
indicating that coupling between bending and 
torsion is very weak and the flutter mode is 
dominated by torsional motion. Therefore, only 
torsion mode is shown for all the calculations. It 
should be noted that these structural properties 
were also used in Refs. 7 and 8. 

The following procedure is used for stability 
prediction: First MISER is run for the selected 
cascade (rotor) using linear theory to identify the 
flutter point. Then MISER is run in conjunction 
with LINFLX2D data base for this flutter condition, 
and the eigenvalues obtained from linear theory 
and LINFLX2D are compared. 

EM Plate Cascade 

To validate the coupling of the unsteady 
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aerodynamic coefficients obtained from LINFLX2D 
with MISER, a nine bladed cascade with a flat plate 
geometry was selected. The flutter results 
obtained are compared with those obtained from 
linear theory. The gap to chord ratio is unity and 
the stagger angle is 45 degrees. For the nine 
bladed cascade, for a given reduced frequency, 
the unsteady force coefficients have to be 
calculated for the possible interblade phase angles 
of 0.0, 40.0, 80.0, 120.0, 160.0, 200.0, 240.0, 
280.0, and 320.0 degrees. 

Using linear theory, for M=0.7 with elastic axis at 
midchord, for the nine bladed cascade mentioned 
above, the instability was observed to occur at k b = 
0.18 and for a phase angle of 80 degrees. 
Therefore, LINFLX2D is run for this value of 
reduced frequency, but for all allowable phase 
angles, and for both plunging and pitching 
oscillations, and a database is prepared. Then 
MISER is run to predict stability boundary with the 
unsteady aerodynamic force coefficient data base 
obtained from LINFLX2D. 

The root locus plot showing the real and imaginary 
parts of the eigenvalues for the flat plate geometry 
obtained for torsion motion is given in Fig. 5 . The 
unsteady force coefficients for eight phase angles 
were obtained with LINFLX2D. For the c =40 
degrees, for which acoustic resonance was 
expected, the LINFLX2D code did not give a 
solution for the grid, oscillation of amplitude, and 
time step considered. This case will be studied 
further in future. It was noted that for a flat plate 
geometry, the computational times were higher 
than those for airfoil geometry (presented below). 
Figure 5 shows the comparison between 
predictions from linear theory and from LINFLX2D. 
It can be seen from Fig. 5 that LINFLX2D compares 
very well with linear theory results for the flat plate 
geometry. The plot shows that the cascade is 
unstable for the properties considered. 

Mistuned Cascade 

Results are presented now for cascades with 
mistuning. The type of mistuning considered is 
the one in which the odd and even numbered 
blades have different torsional frequencies, also 
known as alternate blade mistuning. This has 
been studied in Refs. 7 and 8 with 1% and 1.5% 
alternate mistuning. It is to be noted that in the 
case of 1% mistuning, the frequency ratioaj a /fo 0 


is 1 .005 for all the even blades, and is 0.995 for all 
the odd blades. The study in "Ref. 7 and 8 
considered a 56 bladed cascade, and showed that 
alternate mistuning has a considerable effect on 
stability. A 12 bladed cascade is considered for 
the present study, to limit the computational time. 
MISER was run first for this case as a tuned 
cascade with linear theory, for M=0.7. It was found 
that the cascade is unstable for 1^ = 0.1 and o = 
60°. This value of kj, = 0.1 is used for all the 
calculation presented below. It should be noted 
here that for this cascade acoustic resonance was 
expected for s = 0° and 30° . However, for CIO 
airfoil geometry unsteady solutions were obtained 
with LNFLX2D for these interblade phase angles, 
unlike for the flat plate geometry mentioned in the 
previous section. 

First, the amount of mistuning required to change 
the root locus plot is investigated i.e. whether the 
effects of 1% and 1.5% alternate as reported in 
Ref. 7 also occurs for this cascade. It should be 
noted here that the results from the present 
MISER code were checked with those given in 
Ref.7 for the 56 bladed cascade before starting 
the current calculations. Figure 6 shows the root 
locus plot with 1%, 1.5%, 5%, 10% and 20% 
alternate mistuning for the flat plate geometry for 
M=0.7 and k„ = 0.18. It can be seen that an 
alternate mistuning of 1-10% did not show any 
effect on the eigenvalues. Only 20% alternate 
mistuning showed some effect on the 
eigenvalues. This is in contrast to the level of 
mistuning used in Ref. 7. However, Ref. 7 has 
considered a 56 bladed cascade in incompressible 
flow for which the cascading effect is much 
stronger. 

However, the large amount of mistuning required 
here to have a significant effect on the 
eigenvalues (shown in Fig. 6) is Tiot surprising. 
Similar values of mistuning have been used in Ref. 
8. That study used an amount of 1 to 20% 
mistuning. An examination of Fig. 11 of Ref. 8 
showed that the effect of alternate mistuning 
decreased with increasing Mach number. In light 
of study of Ref. 8 and the present results, it is 
concluded that for the Mach numbers considered 
here a mistuning of more than 20% has to be used 
to show a significant effect of mistuning. 
Therefore the rest of the calculations are done with 
20% alternate mistuning. 

Figure 7 shows the root locus plot obtained for the 
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CIO airfoil for both tuned and 20% alternate 
mistuning, for M=0.7 and k*, = 0.1. Comparison 
with flat plate geometry is also shown. The 
following are noted form Fig. 7: (1) For the tuned 
cascade both flat plate geometry cascade and CIO 
airfoil cascade are unstable for o=60°. However, 
the CIO cascade is less unstable than the flat plate 
cascade. This indicates that for the tuned cascade 
the airfoil steady loading effects are stabilizing, 
compared to that for a flat plate airfoil; (2) the 
addition of 20% alternate mistuning moved the 
cascade with CIO geometry from an unstable to a 
stable condition, where as the flat plate cascade 
remained unstable; (3) It can also be seen that from 
the root locus plot for the CIO airfoil with 
mistuning, two frequencies are identified with 
ct= 240° mode, whereas o=60° mode is completely 
missing. The reason for this behavior is under 
investigation. 

The root locus plot for M=0.8 and k 5 = 0.1 is shown 
in Fig. 8. The following are noted form Fig. 8: 

(1)the tuned flat plate geometry cascade was 
unstable for o=90° and CIO airfoil cascade was 
stable. This indicates that the airfoil steady loading 
effects are stabilizing for the tuned cascade, 
compared to that for a flat plate airfoil; (2) the 
addition of alternate mistuning made the CIO 
geometry more stable, where as the flat plate 
cascade remained unstable; (3) For the CIO airfoil 
it was noted that for tuned case both o=90° and 
c=120° modes have the same amount of damping, 
but with mistuning, a=90° mode became more 
stable than the a=120° mode. 

CONCLUDING REMARKS 

A two dimensional unsteady linearized 
aerodynamic Euler solver, LINFLX2D was used to 
calculate unsteady aerodynamic forces on 
oscillating cascades. The unsteady aerodynamic 
forces were used in the aeroelastic stability code 
MISER to calculate flutter characteristics. Results 
were presented for cascades in subsonic and 
transonic flow. The linearized Euler solver showed 
good correlation with the published data. The 
following were observed during the study. 

(1) the airfoil steady loading stabilized the cascade 
for the Mach numbers considered. 

(2) A relatively high amount of alternate mistuning 
is required to have an effect on the stability of the 
cascade for the Mach numbers considered. 


(3 An accurate steady solution called “low loss 
solution " is required for the unsteady solver 
LINFLX2D to work. This requires a considerable 
amount of computational time. This may be the 
biggest drawback of using this code in a design 
cycle. 

(4) The success of getting an unsteady solution 
from LINFLX2D also depended on the grids used. 

(5) For predicting flutter boundary, calculations 
have to be done for a range of reduced 
frequencies. It is suggested that LINFLX2D data 
base be prepared for three or four reduced 
frequencies, and then use interpolation for flutter 
calculations. These reduced frequency values can 
be selected by first running MISER with linear 
theory. 

(6) The unsteady solution times varied from 20 
minutes to 75 minutes on an SGI workstation. This 
is directly related to the number of iterations 
required for convergence. A solution for flat plate 
geometry took more time than for the CIO airfoil 
geometry. 
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Table 1 


Unsteady Aerodynamic Force Coefficients 

s/c=1.0, stagger=45°, c=180°, k„ =0.5, x0=0.5 
alrf ol 1=1 Oth standard configuration 





plunging 

pitching 


Cl 

a 

-3.56752,1.29455 

-4.244, .73506 


b 

-3.5133,1.2338 

-4.1786, .6609 

flat plate 
M=0.7 


a 

0.55951, -.67647 

.71 385, -.87771 

cm 

b 

0.5297,-0.6417 

.6780, -.8324 



cl 

a 

-3.5675,1.29455 

-4.2439, .73506 


b 

-2.5502,-0.0132 

-2.2566,-. 1537 

airfoil 

M=0.7 


c 

-2.6370,0.0553 

-2.3040, -.0487 


d 

-2.6576,0.0310 

-2.4415, -.0362 


cm 

a 

0.5595, -.67647 

0.71 385, -.87771 



b 

.6497, -.3330 

0.6001,-0.4540 



c 

.61 24, -.3340 

0.5513,-0.4406 



d 

.6346, -.3474 

0.5917,-0.5108 


cl 

a 

-.47916,2.64057 

-.4467,3.10392 


b 

* 

* 

airfoil 

M=0.8 


c 

-2.5945,1.0469 

-1.8953,0.8885 


d 

-2.7240,1.1272 

-2.0532,0.9119 



a 

-.46335,0.03245 

-.78786, -.16579 


cm 

b 

* 

* 



c 

0.3352,-0.6922 

0.2215,-0.6712 



d 

0.3615,-0.6567 

0.2657,-0.6940 


a-linear theory, Ref. 3, b=LINFLX2D (155x41 grid),c=LINFLX2D (141x41 grid) 
d=nonlinear unsteady Euler code, Ref. 17 


could not get solution with this grid 
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Figure 1: Typical section blade model with two degrees-of-freedom. 



Figure 2: Cascade geometry showing stagger angle (6), chord length (c), incidence angle (/) and 
spacing (s ), periodic and injection boundaries used in calculations with multiple grid blocks. 
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distance along the chord, x/c 


Figure 3a: Unsteady pressure difference distribution for a fiat pfate geometry cascade; pitching about 
mid-chord; , s/c = 1.0, 6 = 45°, = 0.7, /= 0°, kb = 0.5, <7=180° 



distance along the chord, x/c 


Figure 3b:Unsteady pressure difference distribution for a subsonic cascade, pitching about midchord' 
s/c =1.0, 61 = 45°, Moo = 0.7, /= 10°.^ = 0.5, <7=180°, a 0 = 2° 
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lower surface 



Figure 4a: Steady Mach number distribution for a transonic cascade; standard configuration No. 10, 
s/c = 1.0, 0 = 45°, Mm = 0.8, / = 1 3°. 



distance along the chord, x/c 


Figure 4b: Unsteady pressure difference distribution for a transonic cascade; standard configuration 
No. 10, pitching about mid-chord; parameters as in Fig. 4a, kb = 0.5, <7=180°, a a = 2°. 
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real part of eigenvalue 

Figure 5: Root locus plot for 9 bladed cascade, torsional motion, flat plate geometry, /cp = 0.1 8, a=0.0, 
p=258.5, r a = 0.5774, structural damping=0.0, M=0.7, /= 10°, cohfco a = 0.357 



Figure 6: Effect of alternate mistuning on eigenvalues in subsonic flow; flat plate geometry, = 0.7, 
/= 10°, /rp = 0.1 ,.N=1 2, pitching about mid-chord; .1%,2%,5%,10%, 20%; structural parameters same 
as in Fig. 5 
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real part of eigenvalue 

Figure 7: Effect of alternate mistuning on eigenvalues in subsonic flow; standard configuration No. 10, 
Mo 0 = 0.7, / = 10°, kfr = 0.5,.N=12, pitching about mid-chord; structural parameters same as in Fig. 5 



real part of eigenvalue 

Figure 8: Effect of alternate mistuning on eigenvalues in transonic flow; standard configuration No. 10, 
Moo = 0.8, / = 13°, kb = 0.1..N = 12, pitching about mid-chord; structural parameters same as in Fig. 5 
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Abst rac t 

A flutter and forced response analysis system based on 
a two-dimensional linearized Euler solver was 
developed. The ASTROP2 code, an aeroelastic 
stability analysis program for turbomachinery, was 
used for this purpose. The ASTROP2 code uses strip 
theory to couple a two dimensional aerodynamic 
model with a three dimensional structural model. 
The code was first modified to include forced response 
capability. The unsteady aerodynamic modeling of 
the ASTROP2 was also changed to accept the loads 
from a linearized Euler solver, LINFLX2D. By 
including the unsteady aerodynamic loads from 
LINFLX2D, it is possible to include the effects of 
transonic flow on flutter and forced response in the 
analysis. The updated code, ASTROP2-LE for 
ASTROP2 code using Linearized Euler aerodynamics, 
is validated by comparing the predictions with those 
obtained using linear unsteady aerodynamic solutions. 

Introduction 

The aeroelastic research program at NASA Glenn 
Research Center is focused on flutter, and forced 
response analysis of turbomachinery. An overview of 
this research was presented in Ref. 1. The review 
showed that a range of aerodynamic and structural 
models have been used to obtain the aeroelastic 
equations. Both time and frequency domain methods 
have been used to obtain unsteady aerodynamic forces 
and to solve the aeroelastic equations. It was noted 
that time domain methods require large computational 
time compared to frequency domain methods, and 
should only be used when non-linearities are expected, 
and when the need justifies the cost. 
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Two approaches were used in obtaining the unsteady 
aerodynamic forces using frequency domain methods. 
In the first approach, Ref. 2, the unsteady 

aerodynamic equations are linearized about a uniform 
steady flow, there by neglecting the effects of airfoil 
shape, angle of attack and thickness. The unsteady 
aerodynamic models developed in Refs. 3-6 based on 
this approach are used in Refs. 7-8 to study the flutter 
and forced response analysis of a compressor rotor 
using a typical section structural model. Some of 
these models were also integrated with a three 
dimensional structural model using strip theory in 
Ref. 9. However, methods developed by this 
approach are restricted to shock-free flows and lightly 
loaded blade rows. 

In the second approach, Ref. 10, the unsteady flow is 
regarded as a small amplitude perturbation about a 
non-uniform steady flow. The unsteady non-linear 
aerodynamic equations are linearized about the non- 
uniform steady flow, resulting in variable coefficient 
linear unsteady aerodynamic equations, which include 
the effects of steady aerodynamic loading due to airfoil 
shape, thickness and angle of attack. Following the 
second approach, Refs. 11-13 developed a steady and 
linear unsteady aerodynamic model based on the 
potential equation. This unsteady aerodynamic model 
was used to study the effect of steady aerodynamic 
loading on flutter stability using a typical section 
structural model in Ref. 14. Subsequently this 
aerodynamic model was integrated with a three 
dimensional structural model using strip theory in 
Ref. 15. 

However, the formulation based on the potential 
equation requires corrections for entropy and flow 
rotation. The Euler equations can be used to correctly 
model rotational and entropy effects associated with 
strong shocks. Unsteady linearized Euler aerodynamic 
models that include the effect of steady aerodynamic 
loading were developed in Refs. 16-18. Recently, a 
two-dimensional linearized Euler code, LINFLX2D, 
was developed under a NASA contract and has been 
reported in reference 19. This code is based on the 
non-linear Euler solver developed in Ref. 20. In Ref. 
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21, the unsteady aerodynamic calculations from 
LINFLX2D were coupled with MISER (Ref. 7); a 
frequency domain aeroelastic stability and response 
code based on a typical section structural model. 
Flutter and forced response calculations were presented 
for cascades in subsonic and transonic flow, with and 
without mistuning. 

The aeroelastic formulation in MISER does not 
represent the behavior of a three dimensional 
structure. An ideal analysis will be to couple a three 
dimensional structural analysis with a three 
dimensional aerodynamic analysis. This may be 
computationally intensive. An intermediate approach 
that couples a two dimensional aerodynamic analysis 
with a three dimensional structural model using strip 
theory, quasi-3D approach, may be less expensive 
computationally. This approach will provide accurate 
results except where three dimensional effects 
dominate. In turbomachines, where the compressors 
and turbines are enclosed, three-dimensional effects 
may not be strong, and strip theory can be used. The 
ASTROP2 code reported in Ref. 9 uses this strip 
analysis approach for aeroelastic stability analysis. 

The primary objective of the present study is to 
develop a quasi-3D aeroelastic code by coupling the 
LINFLX2D code with a three dimensional structural 
model using strip theory. This effort uses the 
ASTROP2 code of Ref. 9 as the basis. The 
ASTROP2 code is a frequency domain stability 
analysis program, and as mentioned before uses strip 
theory to couple a two dimensional aerodynamic 
model with a three dimensional structural model. The 
present version of ASTROP2 has provision for flutter 
calculations, but not forced response. During the 
course of the research effort, the ASTROP2 code is 
extended to include forced response calculation. The 
resulting code is named ASTROP2-LE for ASTROF2 
code using Linearized Euler aerodynamics. Brief 
descriptions of the formulation and method of 
analysis are given in the next section, followed by 
results and discussion. 

Formulation 

The aeroelastic formulation using the linearized 
approach requires solutions from two aerodynamic 
codes and a solution from a structural dynamic 
analysis code. The steady aerodynamic loads are 
obtained from a non-linear Euler solver, NPHASE, 
and the unsteady loads are obtained from LINPLX2D! 
The structural dynamic solution can be obtained from 
any finite element code, analytical solutions or 
measured data. The salient features of both the 
aerodynamic codes, and the aeroelastic formulation are 
described below. 


Aerodynamic model 

Moihlinear Steady Euler Solver. NPHfl Sp. 

The steady aerodynamic model is based on the 
unsteady, two dimensional, non-linear Euler 
equations. The equations in conservative differential 
form are solved in a time dependent body-fitted 
curvilinear reference frame. This transformation 
process and the ensuing numerical method are 
presented in detail in Ref. 20. The equations are 
discretized and solved using a finite volume method 
using a combination of flux difference splitting and 
flux vector splitting scheme. In addition, limiters are 
used to control dispersive errors commonly 
encountered with higher order schemes. The steady 
solutions presented herein are obtained using the basic 
implicit scheme developed in Ref. 20, which is third 
order accurate spatially and second order accurate in 
time. 

Linear Unsteady Euler Solver. LTNFT X2 D 

To obtain the linearized unsteady Euler equations, the 
independent variables in the unsteady non-linear Euler 
equations are expanded in an asymptotic series of the 
form 

U = U(x) + u(x(x,t),t) + higher order terms ( 1 ) 

Assuming the unsteady excitations are harmonic in 
time, and with the first order variable to be 
represented as complex valued, the above equation can 
be written as 

U = U(x) + Re[u(x)exp(i &)t)j + higher order terms (2) 

Here, the term U(x) is of order one and the second 
term is of the order £ . 

Substituting the expansion of Eq. 2 in the non-linear 
unsteady Euler equations, and equating terms of like 
power in £, and neglecting terms of second and 
higher order in £, nonlinear steady equations and 
linear variable coefficient unsteady equations are 
obtained. 

For harmonic blade motions with constant phase 
angle between adjacent blades (interblade phase angle), 
the values of interblade phase angle ( a ) that can 
occur are given in Ref. 22 as 

<J=2 Ttr/N r=0,l,2, N-l (3) 

where N is the number of blades in the cascade. In a 
time domain approach, the number of blocks required 
depends on the interblade phase angle, and small 
phase angles may require large number of blocks to 
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calculate the unsteady aerodynamic forces. However, 
with the linearized approach, the periodic conditions 
are applied on a single extended blade passage region 
i.e., a region of angular pitch, 

6=271 IN (4) 

In solving the linear unsteady equations, the 
independent variables are regarded as pseudo time 
dependent. This allows solutions to be determined 
using conventional time marching algorithms to 
converge the steady and the complex amplitudes of 
the unsteady conservation variables to their steady 
state values. For more details, see reference 19. 

Aeroelastic Model 

As mentioned before, the primary objective of the 
present study is to couple the LINFLX2D solver with 
a three dimensional structural model using strip 
theory. For this effort, the ASTROP2 code of Ref. 9 
is selected as the basis. 

The ASTROP2 code uses strip theory to couple a two 
dimensional aerodynamic model with a three 
dimensional structural model. The present version of 
ASTROP2 can solve only for flutter and not for 
forced response. For clarity and completeness, the 
aeroelastic formulation in ASTROP2 is given here 
with extension to include forced response calculation. 


(l+2i where g is structural damping ratio. The 
expressions for [Aj and {AZ)J are given below. 


In ASTROP2 the blades are divided into strips where 
the aerodynamic forces are calculated. Each strip has 
two degrees of freedom, a plunging displacement, h, 
motion perpendicular to chord, and a pitching 
(torsion) displacement, a, rotation about the leading 
edge of the strip (Fig. lb). 

Using the normal modal values obtained from a free 
vibration analysis, the equivalent h and a for k lh strip 
of the s‘ h blade are given as summation of normal 
modes as 



OC x CL 2 *.-‘6 L nm 




As mentioned before, an N blade cascade can vibrate 
in N interblade phase angles given by Eq. 3. The 
aerodynamic forces at each strip for the V th interblade 
phase angle mode can be written as a sum of motion 
dependent and motion independent forces as 

[ L rk] = [L\r k ] + {L2 rk ) (7.1) 

where 


ASTROP2 uses the normal mode approach for 
aeroelastic analysis. The equations of motion for the 
s* blade of the cascade, Fig. la, can be written as 


[ii„i= 


^hhr ^har 
Jodir ^attr 


M - 


Khr har 
. ^cxhr l acxr J 




K*]fe} + K, ]{<?,} = Klfe} + {ADj (5) 

where and [ K gx j are generalized mass and 

stiffness matrices, which are diagonal, jgj is the 
generalized displacement vector, and is the 

motion dependent generalized aerodynamic load 
matrix, and {AZ^} is the motion independent 

generalized aerodynamic load vector. The motion 
dependent forces cause flutter, and motion independent 
forces cause forced response (forced vibration). The 

matrices [a^], and [Aj are of NM x NM 

size, {*,} and { AD s} are of NM x 1 size, where 
NM is the number of normal modes used in the 
analysis. The matrices and £ K^ s j are given 

by a free vibration analysis. It should be noted that 
in the stiffness matrix, the structural damping could 
be added by multiplying the stiffness coefficients by 


{ L2 rJ = jr r } (7-2) 

[war J 

where the subscript ‘w* refers to force from wake i.e. 
motion independent aerodynamic forces. The values 

° f hhr' har' ^ ahr ' har’ Lhr ' 30(1 Lor for a £ iven 

interblade phase angle are calculated from the unsteady 
aerodynamic forces obtained from linear theory, Ref. 
4 or LINFLX2D, Ref. 19. 

The generalized aerodynamic matrices [A vr ] and 
|AZ) ?r | for the entire blade are given by 

K]=jK] 7 ’[n rit ]^ = JK,K (8.i) 

o 

{*»„} = (8.2) 

0 

where L is the blade length and T indicates transpose. 
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Integrals in Eq. 8 are evaluated using numerical 
integration as 

nstrip 

Kr] = XKrtfe (9.1) 

k=\ 



For a stability calculation (flutter), the motion- 
independent forces (AD) are set to zero and the 
eigenvalue problem is obtained as 


NSTRIP 

{ AZ) "}= (9.2) 

k=l 

where ds k is the strip width, and NSTRIP is the 
number of strips the blade is divided into. 

For a tuned cascade, in which all the blades are 
identical structural properties, the interblade phase 
angle modes are uncoupled. In this case, the 
aeroelastic modes consist of individual blades 
vibrating with equal amplitudes with a fixed 
interblade phase angle between adjacent blades. 
Hence, for this problem, the motion of the typical 
blade is written as 

k}={?o,F“ do) 


[M-rfe]]fc,Mo} ( i5) 

The solution of the above eigenvalue problem results 
in NM complex eigenvalues of the form 


i ~ 

The real part of the eigenvalue (fl ) represents the 
damping ratio, and the imaginary part (V ) represents 
the damped frequency; flutter occurs if ]± > 0 for any 
of the eigenvalues. The eigenvalue problem is solved 
for each of the N permissible values. The critical 
phase angle is identified as the one that results in the 
lowest flutter speed. 


j=*'Vr =n+iv 


(16) 


The equations of motion for the blade for each Th f aeroeIastic res P onse of the blades induced by 

interblade phase angle can be written as Wakes 1S calculated fr °m equation ( 1 3) as 


= °> 2 K] k r y ( “ +<V) +CQ 1 [AD r ]e i( “ +a ' s) 

(11) 

The matrices KJ’ Kir]- K] are of NM x 

NM size, KJand {&D r } are of NM x 1 size. 

Since the blades are identical, the same equation is 
obtained for each blade; equation 1 1 can be solved for 
N different values of the interblade phase angle given 
by Eq. 3. 

Dividing both sides by an assumed frequency,^, 

dropping the subscripts, and rearranging, the 
equations can be written as 

Vi P ]-y[Q][{qa r } = Y{AD} (13) 

where 


CO Q 

[Q] = [[M} + [A]] 


{*«■}= -[[ p ]-r[e]r'r{AZ9} a 7) 

Stability calculation procedure 

The aerodynamic coefficients are calculated before the 
eigenvalue problem can be set up and solved. Since 
the unsteady aerodynamic coefficients depend on the 
frequency of oscillation, it is necessary to assume a 
frequency CD 0 (actual input to the code may be the 
reduced frequency of blade vibration based on chord, 
K = COC / U , where U is the free stream velocity) in 
advance to be able to calculate the aerodynamic 
coefficients. The aerodynamic coefficients are 
functions of inlet Mach number M y and interblade 
phase angle (J r , in addition to cascade geometric 
parameters. In the present study, for a given inlet 
Mach number, the reduced frequency is varied until 
the real part of one of the eigenvalues fl becomes 
zero while the real parts of the remaining eigenvalues 
are negative. The assumed flutter-reduced frequency 
k c j and the calculated flutter frequency V f are both 

based on CQj . Thus, these two can be combined to 

eliminate Q)j and the flutter speed is obtained, 

namely, V jC0) Q lk L ^ . Since the inlet Mach number 

is known, this flutter speed gives the inlet condition 
(speed of sound, aoo ) at which the cascade will be 
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neutrally stable for the given Mach number. This 
procedure can be repeated to obtain a plot of flutter 
speed versus Mach number. Knowing the operating 
conditions, it is possible to determine whether flutter 
will occur within the operating region and if so, the 
Mach number and frequency at flutter. 

Analysis procedure 

The flutter and forced response analysis in 
ASTROP2-LE consists of running five codes, (I) a 
structural dynamic analysis code, (II) 2DSTRIP (El) 
NPHASE, (IV) LINFLX2D, and (V) 2DASTROP. It 
is to be noted that ASTROP2 code is a combination 
of 2DSTRIP and 2DASTROP codes. The analysis 
procedure is explained in five steps. In step 1 a 
vibration analysis is performed for the blade. The 
output is natural frequencies and mode shapes. This 
output is used by 2DSTRIP. In step 2, strips are 
selected, and 2DSTRIP is run to calculate relative 
Mach numbers, sweep angles, stagger angles, chord 
values, and strip widths at these strips. During this 
run, the three dimensional modal values are also 
interpolated at each strip, and equivalent pitching and 
plunging modal values are calculated. In step 3, a 
steady aerodynamic solution at these strips is obtained 
from NPHASE. The steady aerodynamic solution is 
written as a database. Step 4 consists of running 
LINFLX2D ^ for assumed number of reduced 
frequencies, interblade phase angles, for pitching and 
plunging modes. The analysis is carried out for unit 
amplitude of vibration for all the strips, and the 
unsteady aerodynamic solutions are stored as a 
database. In step 5, the 2DASTROP code is executed 
to calculate flutter using the eigenvalue approach and 
to calculate forced response. 

St rategy, to reduce the, number o f calculatinnc 

The main purpose of the present effort is to accurately 
couple the unsteady aerodynamic solutions from 
LINFLX2D to ASTROP2 code. In order to reduce 
the number of calculations, the following 
observations are made on the number of calculation 
required for flutter. 

Let each blade be divided into NSTRIP number of 
strips. In general, for a given inlet Mach number, the 
flow conditions, (relative Mach number and angle of 
attack), and the geometric conditions, (gap to chord 
ratio, stagger angle, and airfoil shape) are different at 
these strips. Therefore, the number of steady 
aerodynamic solutions (number of NPHASE runs) 
required is (NSTDY = NSTRIP * NMACH ), where 
NMACH is the number of inlet Mach numbers to be 
considered in the study. By considering a straight, 
untapered stator blade, the flow and geometric 
conditions will be same at all strips, reducing the 


number of steady aerodynamic calculations to 
NMACH, i.e. NSTD Y=NM ACH . 

In the present study two inlet Mach numbers will be 
considered i.e. NMACH=2. Therefore the number of 
study runs is NSTDY = 2*NSTRIPS. Assuming 
identical flow and geometric conditions, the number 
of NPHASE (study solution) runs reduces to NSTDY 
= 2 . 

In the case of the unsteady aerodynamic solution, two 
other parameters, reduced frequency ( k c ) and 
interblade phase angle ( a ) have to be considered. If 
NREDF is the number of reduced frequencies, and 
NSIGMA is the number of interblade phase angles, 
then the number of unsteady solutions for flutter 
(number of LINFLX2D runs) is, NUSTDY = 
NMACH *NSTRIP *NSIGMA * NREDF. However, 
if same geometric and flow conditions are assumed at 
all strips, NUSTDY is reduced to NUSTDY = 
NMACH*NSIGMA*NREDF. It is to be noted here 
that NSIGMA is equal to the number of blades of the 
cascade. 

For NM ACH=2, the number of unsteady runs 
NUSTDY= 2*NSTRIP*NREDF*NSIGMA. Again 
assuming identical flow and geometric conditions, 
NUSTDY=2*NREDF*NSIGMA. 

Results a nd Discussion 

Results presented here are meant to show that the 
analysis procedure given in the previous sections has 
been implemented correctly in LINFLX2D and 
ASTROP2. Therefore, calculations are made for a 
non-rotating cantilevered blade, representing a stator 
blade of a turbomachinery component. Results are 
presented for an assembly of 12 tuned blades on a 
rigid disk coupled only aerodynamically. Similar 
geometry was considered in Ref. 23. Two inlet Mach 
numbers are considered M=0.7 and 0.8. Two airfoils, 
a flat plate and the tenth standard configuration, CIO 
airfoil cross sections. Ref. 25, are considered. The 
gap to chord ratio (s/c) is 1.0 and the stagger angle is 
45 degrees. For the CIO airfoil and for M=0.8, the 
flow is transonic with a shock near the quarter chord. 

Code Validation: 

ASTROP2 code: 

The ASTROP2 code was validated for flutter 
prediction in Ref. 9 by calculating the flutter 
boundary of the SR3CX2 advanced propeller. The 
vibration characteristics of the propeller blade were 
obtained from the finite element structural analysis 
code NASATRAN. Linear theory of Ref. 5 was used 
for calculating the unsteady aerodynamic forces. The 
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predicted flutter Mach number was conservative 
compared to the measured value. For more details on 
the ASTROP2 code, see references 9, 24 and 25. 

LINFLX2D code: 

To validate the LINFLX2D code, unsteady 
aerodynamic pressure differences were calculated for a 
flat plate cascade, and a cascade comprising of CIO 
airfoils with inflow Mach numbers of M=0.7 and 
M— 0.8, for a reduced frequency of oscillation based on 
chord (k c ) of 1.0. The unsteady pressure difference 
is non-dimensionalized by (upper surface pressure- 
lower surface pressure) / (airdensity*U**2*ampIitude 
of oscillation ). An H-grid of 141x41 grid is used for 
the study. There are 80 points on the airfoil, and the 
inlet boundary is 5 chords from the airfoil leading 
edge and the exit boundary is 10 chords from the 
airfoil trailing edge. Steady aerodynamic solutions 
were obtained for these Mach numbers usin< J 
NPHASE before running LINFLX2D. 

Subsonic Inflow 

Figure 2a shows the steady Mach number distribution 
obtained for the CIO airfoil cascade for an inlet Mach 
number, M=0.7. The flow to the cascade is at 10 
degrees angle of attack. It can be seen that for this 
Mach number the flow is shock free. The unsteady 
pressure difference distribution for an interblade phase 
angle, (7, of 180 degrees is shown in Fig. 2b. The 
blades are oscillating in pitch about the midchord. 
Figure 2b shows predictions from linear theory ( Ref. 

4 ) for flat plate geometry, from the nonlinear Euler 
(Ref. 20) and the linearized Euler code, LINFLX2D 
for the CIO airfoil geometry. The predictions from 
L1NFLX2D correlate well with the nonlinear Euler 
results. As expected, the results show significant 
difference with linear theory, because the effects of 
airfoil geometry, and angle of attack are not included 
in the linear theory. 

Figure 2b also shows the unsteady pressure difference 
distribution obtained with LINFLX2D for the flat 
plate geometry, and its correlation with linear theory. 

It can be seen that for the flat plate geometry results 
from linear theory and LINFLX2D match very well. 

At M=0.7 the flow field is linear; hence, linear theory 
and the linearized Euler results are expected to 
correlate well. This validates the LINFLX2D code for 
subsonic flows. 

Transonic Inflpw 

For the CIO airfoil cascade geometry, and for an inlet 
Mach number of 0.8, the flow is transonic with a 
normal shock occurring in each blade passage. The 
steady Mach number distribution obtained from 


NPHASE is shown in Fig. 3a for a steady angle of 
attack of 13 degrees. A normal shock occurs on the 
suction surface at about 28% of the chord. Results 
from the full potential solver, Ref. 12, are also 
included for comparison. Both results agree well, 
except that the nonlinear steady Euler solver' 
NPHASE, predicts the shock location slightly 
downstream of that predicted by the full potential 
solver of Ref. 12. A shock fitting procedure was 
used in Ref. 12, whereas the shock is captured 
naturally in NPHASE. 

The unsteady pressure distribution is shown in 
Fig. 3b, along with a comparison with linear theory 
(Ref. 4) and the unsteady nonlinear Euler solver. The 
unsteady results are for pitching about mid-chord with 
~ TO. O’ = 180° and amplitude of oscillation, 
OC 0 , = 2°. As expected, linear theory does not show 

any shock, indicating that unsteady analysis based on 
linear theory will not be accurate for cascade flutter 
analysis in transonic flow. Some discrepancy is seen 
between the unsteady results obtained from 
LINFLX2D and from the unsteady nonlinear Euler 
solver. For the transonic flow case the following 
factors may have more effect than for in subsonic 
case: (1) grids used, (2) accuracy lost when time series 
data is processed to obtain frequency data or (3) the 
steady solution on which LINFLX2D based is not a 
low loss solution. Even though these issues need 
separate study. Fig. 3b shows a fair correlation 
between LIFLX2D results and non-linear unsteady 
Euler. 

Similar results were obtained for plunging motion for 
both subsonic and transonic inflows. 

Stability and Res ponse Calculations 

As mentioned before, a 12 blade stator is considered 
for stability and response calculations. Since the 
main aim of this paper is to demonstrate the 
ASTROP2 and LINFLX2D code coupling, analysis is 
earned out for an assumed geometry of straight, 
untwisted stator blades at two inlet Mach numbers, 
0.7 and 0.8. This reduces the number of steady runs, 
NSTDY, 2 (one for M=0.7 and one for M=0.8) and 
the number of unsteady runs, NUSTDY, to 
24*NREDF for a 12 blade cascade. The same grids 
that were used for the unsteady validation and the 
NPHASE steady solution obtained in the previous 
section were used again. 

The blade is divided into 10 strips (NSTRIP=10). 
Table 1 shows the aerodynamic input parameters, and 
table 2 shows the modal values used in this study. 
These are same as those used in Ref. 23. The modal 
values at the strips are used for generalized force 
calculation. Two modes are used in the analysis. 
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The First mode is a pure torsion mode and the second 
mode is a coupled bending-torsion mode. The natural 
vibration frequencies are 81.376HZ and 148.02 HZ 
respectively. 

Stability calculating 

The unsteady aerodynamic force coefficients are 
calculated using LINFLX2D for harmonic blade 
vibration in plunging and pitching modes for a 
reduced frequency, k c , of 0.2 based on chord. The 
calculated unsteady force coefficients are used by 
2DASTROP code to calculate the elements of [A], 
With identical flow and geometric conditions, the 
only contribution to the generalized aerodynamic force 
matrix is due to the different modal values at different 
strips. The reduced frequency used to calculate the 
unsteady aerodynamic force coefficients, k c =0.2, 
corresponds to a frequency of 96.07 HZ. 

Figure 4 compares the root locus plot of the 
eigenvalues obtained from LINFLX2D with that 
obtained using linear theory unsteady aerodynamics 
for a flat plate geometry. The plot is shown for the 
first mode for M=0.7. Good correlation is observed. 
At M=0.7 the flow field is linear hence, linear theory 
and linearized Euler are expected to correlate well with 
each other. The linear theory unsteady aerodynamic 
code was an integral part of ASTRO P2, and 
LINFLX2D is the new code that is coupled with 
ASTROP2. The root locus plot correlation shows 
that the coupling of the L3NFLX2D unsteady 
aerodynamic database to ASTROP2 is accurate. It is 
to be noted that since the first mode is a pure 
torsional mode, the frequency of the aeroelastic 
system is close to the first natural frequency of the 
blade. 

To explore the effects of airfoil shape and transonic 
flow the analysis was carried out with the CIO 
geometry and for M=0.7 and 0.8. Figure 5 shows the 
root locus plot for the first mode at an inflow M of 
0.7 for the flat plate and CIO airfoil geometries. It 
can be seen that the blade is more unstable when the 
effect of airfoil geometry is included. 

The root locus plot for the first mode for M=0.8 is 
shown in Fig. 6. Here the root locus plot shows the 
effect of both airfoil shape and transonic flow effects. 

It is seen that the blade is slightly more unstable than 
in figure 5 when the effects of airfoil geometry, and 
shock are included. 


The response of the blade to a vortical disturbance is 
calculated for flat plate geometry and CIO geometry 
and compared with that obtained from linear theory. 


In general the response calculation can be done by 
including the motion dependent aerodynamic matrix, 
[A], which adds the contribution of aerodynamic 
damping to the response. The forced response is 
usually calculated for an aeroelastically stable system 
(positive aerodynamic damping), and the contribution 
of [A] to response can be ignored. For the 
calculations presented here, the contribution of [A] is 
neglected. However, a structural damping ratio of 

0. 002 is added. 

Figure 7 shows the tuned aeroelastic response for r=6 

1. e. o =180 degrees for flat plate geometry. The 
unsteady aerodynamic coefficients are obtained at 
kc = 1*0 for M=0.7. Calculations are done with a 
155x41 grid for which the inlet and exit boundaries 
are located at one chord from the leading and trailing 
edges respectively. There are 55 points on the airfoil 
for this grid. The forcing frequency range investigated 
is limited to a small range around the 1 st mode 
frequency. For the tuned cascade all the blades will 
have equal amplitudes. The amplitude of response is 
a function of the frequency ratio, oo / co t) . Figure 7 
shows the 1st generalized degree of freedom (ql) 
response obtained using linear theory and present 
LINFLX2D code. It can be seen that calculations 
from linear theory and from LINFLX2D are identical 
indicating that the coupling of LENFLX2D 
coefficients to ASTROP2 code for response 
calculations is also accurate. 

Figure 8 shows the ql response obtained for the CIO 
airfoil for M=0.7. Comparing the response of CIO 
with that of flat plate, it can be seen that the response 
has increased about 140%. This is due to high steady 
loading on the airfoil. 

Figure 9 shows the ql response obtained for the CIO 
airfoil for M=0.8. Comparing the response of CIO 
with that of flat plate, it can be seen that the response 
has increased to only about 68%. Observing from 
Fig 8 that the airfoil shape increased the response to 
about 140%, the decrease in the response in Fig 9 can 
be attributed to the presence of shock. However, the 
flow being transonic, it is planned to check the 
present calculations by running LINFLX2D for a 
different grid. 

Concluding Remark.*? 

The unsteady linearized Euler aerodynamic solver, 
LINFLX2D, has been successfully coupled with the 
ASTROP2 aeroelastic analysis code. The resulting 
code, ASTROP2-LE, is validated by comparing 
predictions from linear theory for flat plate geometry. 
Comparison was done for both flutter and forced 
response. Results were also presented for cascades in 
subsonic and transonic flow for a standard cascade 
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section known as CIO geometry. The following were 
observed during the study. 

(1) The steady loading due to the airfoil shape, angle 
of attack destabilized the cascade for the Mach 
numbers considered. 

(2) The steady loading due to the airfoil shape, angle 
of attack increased the blade response for subsonic 
flow and decreased response for transonic flow 
compared to that for a flat plate. 

(3) An accurate steady solution called “low loss 
solution “ is required for the unsteady solver 
LINFLX2D to work. This requires a considerable 
amount of computational time. This may be the 
biggest drawback of using this code in a design cycle. 

(4) The success of getting an unsteady solution from 
LINFLX2D also depended on the grids used. 

(5) It is to be noted that the number of LINFLX2D 
solutions required is directly related to number of 
sfrips. Care has to be taken in selecting the number 
of strips to reduce the number of calculations. 

(6) For predicting the flutter boundary, calculadons 
have to be done for a range of reduced frequencies It 
is suggested that a LINFLX2D database be prepared 
for three or four reduced frequencies, and then 
interpolation be used for flutter calculations. These 
reduced frequency values can be selected by first 
running ASTROP2-LE with linear theory. 

(6) The unsteady aerodynamic solution times varied 
from 20 minutes to 75 minutes on an SGI 
workstation. This is directly related to the number of 
iterations required for convergence. A solution for 
flat plate geometry took more time than for the CIO 
airfoil geometry. 

(7) The LINFLX2D code can be used with any 
aeroelastic code that uses a typical section aeroelastic 
model as its basis. 
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Tablel: Aerodynamic Input Parameters at the Strips 

Atmospheric pressure (psi)=l 3. 1023; Speed of sound (fps)=l 130.0 


Strip Index 

Stagger Angle 
(Degrees) 

Chord Length 
(Inches) 

Gap/Chord 

Ratio 

Radius 

(Inches) 

Strip width 

1 

45.0 

3.145 


5.407 

0 1 330 

2 

45.0 

3.145 

1.0 

5 540 

n 1 33^ 

3 

45.0 

P 3.145 

1.0 

5.674 

kJ • 1 J J J 

0 1335 

4 

45.0 

3.145 


5.807 

0 1335 

5 

45.0 

3.145 

1.0 

5.941 

0 1335 

6 

45.0 

3.145 

1.0 

6.074 

v/. Jl J J J 

0 1 335 

7 

45.0 

3.145 


6.208 

0.1335 

8 

45.0 

3.145 


6.341 

0 1335 

9 

45.0 

3.145 1 

1.0 

6.475 

v. 1 J J J 

0 1 335 

10 

45.0 

nn i 

1.0 

6.608 

V. I J J J 

n n^o 




u.uooy 
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Table 2: Mode shapes and frequencies used 
in the study 

Table 11 : Mode 1: natural frequency -8 1376 HZ 


Table 12: Mode 2: natural frequency^ 148.02 HZ 


Strip Index 


Torsion 

i 

0.0 

0.328 


0.0 

0.444 

3 

0.0 

0.584 

4 

0.0 

0.627 

5 

0.0 


6 

0.0 

0.742 

7 

0.0 

mi —i 

8 

0.0 

0.889 

9 

0.0 

0.937 

10 


1.000 


Strip Index 

— 

Torsion 

1 

-0.082 

0.329 

2 

-0.087 

0.446 

3 

-0.070 

0.586 

4 

-0.058 

0.630 

5 

-0.034 

0.696 

6 

0.0 

0.745 

7 

0.076 

0.796 

8 

0.203 

0.892 

9 

0.358 

0.940 

10 

0.502 

1.003 



Figure la. : ASTROP2 coordinate system for a 
rotating blade 



Figure lb. : Section A-A showing rigid pitching (a) 
and plunging (h) motions for the strip ( reference axis 
=Ieading edge) 



















































Fig2a.Steady Mach number distribution for a 
subsonic cascade, CIO airfoil, s/c=1.0, staggei=45 
degrees, M=0.7, angle of attack = 10 degrees. 



Figure 3a. Steady Mach number distribution for a 
transonic cascade; CIO, sic = 1.0, stagger =45° 

M oo = 0.8, i = 13°. 



Figure 2b: Unsteady pressure difference distribution 
for a subsonic cascade, pitching about midchord, 
s/c= 1 .0,stagger=45 deg., M=0.7, kc=1.0,sigma=180 
deg., alfaO=2 deg., 1=0 for flat plate and 1=10 for CIO 
airfoil 



Figure 3b: Unsteady pressure difference distribution 
for a transonic cascade; CIO configuration, pitching 
about mid-chord ° 
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Figure 4. Root locus plot for 12 blade cascade, first mode, flat plate geometry 

kb=0. 1, structural damping=0.0,M=0.7 
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Figure 5. Root locus plot for 12 blade cascade, first mode, flat plate and 
CIO geometries, kb=0. 1, structural damping=0.0,M=0.7 
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real part of eigenvalue 


Figure 6. Root locus plot for 12 blade cascade, first mode, flat plate and 
CIO geometries, kb=0.1, structural damping=0.0,M=0.8 



Figure 7. First mode response for a tuned cascade, flat plate geometry 
structural damping =0.002, kb=0.5,N=12, flat plate geometry ,M=0.7’ 
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Figure 8. First mode response for a tuned cascade, flat plate geometry 
structural damping =0.002, kb=0.5,N=12, flat plate and CIO geometries, M=0.7 



Figure 9. First mode response for a tuned cascade, flat plate geometry, 
structural damping =0.002, kb=0.5,N=12, flat plate and CIO geometries,M=0.8 
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Ab stra ct 

This paper describes the development of LINFLUX- 
AE, a turbomachinery aeroelastic code based on the 
linearized unsteady 3D Euler solver, LINFLUX. The 
steady solution, required by LINFLUX is obtained 
from the nonlinear Euler/Navier Stokes solver 
TURBO-AE. The paper briefly describes the salient 
features of TURBO-AE, LINFLUX and the details of 
the aeroelastic extension. The aeroelastic 
formulation is based on a modal approach. The 
unsteady aerodynamic forces required for flutter and 
forced response are obtained by running LINFLUX 
for each mode for each interblade phase angle and for 
a given frequency. An eigenvalue formulation is 
used for flutter analysis. The forced response is 
calculated from the modal summation of the 
generalized displacements. The NASA / GE Energy 
Efficient Engine (E- cubed) fan configuration is 
chosen as a test case and the results are compared 
with those obtained from TURBO-AE and 
ASTROP2. 

Introduction 

An overview of the aeroelastic analysis methods for 
turbomachines, Ref. 1, shows both time and 
frequency domain methods have been used to obtain 
unsteady aerodynamic forces and to solve the 
aeroelastic equations. It was noted that time domain 
methods require large computational times compared 
to frequency domain methods, and should only be 
used when non-linearity’s are expected, and when the 
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need justifies the cost. Two approaches were used in 
obtaining the unsteady aerodynamic forces using 
frequency domain methods. In the first approach, 
Ref.2, the unsteady aerodynamic equations are 
linearized about a uniform steady flow, thereby 
neglecting the effects of airfoil shape, incidence and 
thickness. The unsteady aerodynamic models 
developed in Refs. 3-6 based on this approach were 
used in Refs. 7-8 to study the flutter and forced 
response analysis of a compressor rotor and propfans. 
However, methods developed by this approach are 
restricted to shock-free flows through lightly loaded 
blade rows. 

In the second approach, Ref. 9, the unsteady non- 
linear aerodynamic equations are linearized about the 
non-uniform steady flow, resulting in variable 
coefficient linear unsteady aerodynamic equations, 
which include the effects of steady aerodynamic 
loading due to airfoil shape, thickness and angle of 
attack. Following the second approach, steady and 
linear unsteady aerodynamic models have been 
developed based on the potential equation , Refs. 10- 
12. This unsteady aerodynamic model was integrated 
with a three dimensional structural model using strip 
theory in reference 13. 

The formulation based on the potential equation 
requires corrections for entropy and rotational flow. 
The Euler equations can be used to correctly model 
rotational and entropy effects associated with strong 
shocks. Unsteady linearized Euler aerodynamic 
models that include the effect of steady aerodynamic 
loading were developed in Refs. 14-16. Recently, 
two and three-dimensional linearized Euler codes 
named respectively LINFLUX2D and LINFLUX 
were developed in references 17 and 18 under a 
NASA contract. These codes were based 
respectively on the non-linear Euler solver developed 
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in Ref. 19 and Refs. 20-21. In Ref. 22, the unsteady 
aerodynamic calculations from LINFLUX2D were 
used with MISER (Ref. 7), an aeroelastic stability 
and response code based on a typical section 
structural model. Flutter and forced response 
calculations were presented for cascades in subsonic 
and transonic flow, and with and without mistuning. 
However, the aeroelastic formulation with MISER 
does not represent the actual behavior of a three 
dimensional structure. An ideal analysis will be to 
couple a three dimensional structural analysis with a 
three dimensional aerodynamic analysis. 

The primary objective of the present study is to 
develop a 3D-aeroe1astic code by coupling the 3-D 
linearized Euler analysis code, LINFLUX of Ref. 1 8, 
with a three-dimensional structural dynamics model. 
A modal approach will be used to do the flutter and 
response analysis. The resulting code will be called 
LINFLUX-AE to indicate the AeroEIastic extension 
to LINFLUX code. Brief descriptions of the 
formulation and method of analysis are given in the 
next section, followed by results and concluding 
remarks. 

For m ul ati on 

The aeroelastic formulation involves solutions from 
three models: a steady aerodynamic model, a 
linearized unsteady aerodynamic model and an 
aeroelastic model coupling the unsteady aerodynamic 
solution with the structural dynamic solution. The 
aerodynamic models and the aeroelastic formulation 
are described in this section. 

Aerodynamic models 

Non-linear Steady Euler Solver 

The steady aerodynamic loading required for the 
linearized solver, LINFLUX, is obtained from the 3D 
Euler / Navier-Stokes solver TURBO-AE, Ref, 23. A 
very brief description of the TURBO code is 
presented in this section. Additional details 
regarding the code are available in Refs. 20-21 and 
23. 

The TURBO code was originally developed in Ref. 
20 as an inviscid flow solver for modeling the flow 
through multistage turbomachinery. It has the 
capability to handle multiple blade rows with even or 
uneven blade count, stationary or rotating blade rows 
and blade rows at an angle of attack. Multiple blade 
passages are included in the calculation, when 


required. Additional developments were made in 
Ref. 21 to incorporate viscous terms into the model. 
The code can now be applied to model realistic 
turbomachinery configurations with flow phenomena 
such as shocks, vortices, separated flow, secondary 
flows, and shock and boundary layer interactions. 

The code is based on a finite volume scheme. Flux 
vector splitting is used to evaluate the flux Jacobians 
on the left hand side of the governing equations and 
Roe's flux difference splitting is used to form a 
higher-order TVD scheme to evaluate the fluxes on 
the right hand side. Newton sub-iterations are used at 
each time step to maintain higher accuracy. A 
Baldwin-Lomax algebraic turbulence model is used 
in the code. Only the Euler option of the TURBO 
code is used in the present study. 

Linear Unsteady Euler Solver. LINFLUX 

To obtain the linearized unsteady Euler equations, 
Ref. 18, the dependent variables in the unsteady non- 
linear Euler equations are expanded in an asymptotic 
series of the form 

U= £/(.*)+ u(x(xj),t)+ higher order terms (1 ) 

Assuming the unsteady excitations are harmonic in 
time, and with the first order variable to be 
represented as complex valued, the above equation 
can be written as 

U= U ( x) + Re [u( x) exp (icot) } +H.O.T. (2) 

Here, the term U(x) is of order one and the second 
term is of the order € . 

Substituting the expansion of Eq. 2 in the non-linear 
unsteady Euler equations, and equating terms of like 
power in £ , and neglecting terms of second and 
higher order in £, nonlinear steady equations and 
linear variable coefficient unsteady equations are 
obtained. 

For harmonic blade motions with constant phase 
angle between adjacent blades (interblade phase 
angle), the values of interblade phase angle (<T) that 
can occur are given as (Ref. 24). 

O r =2ltr! N ; r = 0,1,2 N-l (3) 

where N is the number of blades in the blade row. In 
a time domain approach with periodic boundary 
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conditions, the number of blade passages required to 
be modeled depends on the interblade phase angle, 
and small phase angles may require large number of 
blade passages to calculate the unsteady aerodynamic 
forces. However, with the linearized approach and 
the assumed harmonic variation in time, the periodic 
conditions are applied on a single extended blade 
passage region i.e., a region of angular pitch, 

6 = 2x/N (4) 

In solving the linearized unsteady equations, the 
dependent variables are regarded as pseudo-time 
dependent. This allows solutions to be determined 
using conventional time -marching algorithms to 
converge the steady and the complex amplitudes of 
the unsteady conservation variables to their steady 
state values. For more details, see reference 1 8. 

Aeroelastic model 

The aeroelastic equations of motion for each blade of 
the blade row of a turbomachine can be written as 

[M]{X} + [*]{X} = (P(X),/)}+{F(f)} (5) 

where [M] is the mass matrix, [K] is the stiffness 
matrix that includes the effects of rotation 
(centrifugal stiffening and softening), {X} is the 
vector of blade deflections at the finite element grid 
points, {P(X),t} is aerodynamic force vector due to 
blade vibration that may cause instability, and (F(t)} 
is the aerodynamic force vector due to gust leading to 
fatigue failures due to forced response. 

A modal formulation is used to proceed further. The 
starting point is to obtain a free vibration solution of 
the blades giving generalized masses, natural 
frequencies, and mode shapes. This solution can be 
obtained either by using commercial finite element 
codes or by writing one’s own analysis solver. The 
general vibratory motion then can be expressed as a 
superposition of the contributions of the normal 
modes, [</>]= [0p0 2 , <(> NM ], as 

i=NM 

{X} = [(t>] T {q) = (6) 

1 = 1 

where {#} is the generalized displacement, and NM 

is the number of normal modes to be used in the 
analysis. Substituting Eq. 6 in Eq. 5, and post 
multiplying the result by [$>] 7 , the equations of 
motion for the s lh blade can be written as 


[ M .]fe} + [^K«.}=[4]{?,} + {AD,} (7.1) 

where 

= (72) 

[4]=M r {W0,(> : 

{AD,} “WW')} 

where [AfJ and \K S ] are generalized mass and 
stiffness matrices, which are diagonal, is the 

generalized displacement vector, and [Aj is the 
motion dependent generalized aerodynamic load 
matrix, and {AX^j is the motion independent 

generalized aerodynamic load vector. The motion 
dependent forces cause flutter, and motion 
independent forces cause forced response (forced 
vibration). The matrices [M,], [£,]. and [/),] are 

of NM x NM size, and [AD^} are of NM x 1 
size. 

For a tuned cascade (or rotor), in which all the blades 
are identical, the aeroelastic modes consist of 
individual blades vibrating with equal amplitudes 
with a fixed interblade phase angle between adjacent 
blades. Hence, the motion of the sth blade in rth 
interblade phase angle mode can be written as 

{?,} = k.K = {«.,}<'" «"■* ® 

where (J r is given by Eq. 3. 

Thus the equations of motion for the blade, Eq. 7.1 
becomes 

[K,] {. q „}e ««•*> 

= [4] {AD,} e*™* (9) 

Since the blades are identical, the same equation is 
obtained for each blade. Thus, no additional 
information can be obtained by assembling the 
equations for all the blades on the disk. Instead, 
equation (9) can be solved for N different values of 
the interblade phase angle given by equation (3). 


Dropping the subscript s, since each blade is 
identical, and canceling out the exponential terms, 
Eq. 9 can be written as 

-® 2 [m] {y} ♦ M{y}-[/t r ]{y}={Az>,}<io> 

where { Y) is same as 

Calculation of Elements of [ArJ and [<p] 

The elements of [A/ 5 ], [A^j and [ <p ] are given by 

the free vibration analysis using a commercial code 
like ANSYS/ NASTRAN or by writing own analysis 
program. The matrices [Af 5 ] and are diagonal. 

The elements of [A/J and [A^] are related as 

Kjj *= M u (Dy 2 ( 11 ) 

where CO <v is the natural frequency of the i lh mode. It 

is to be noted that usually the mode shapes are mass 
normaized giving M H = 1.0. 

Calculation of elements of [A r ] and {AD r j 

The elements of [A r ] and |AZ) r | are obtained from 

LINFLUX code. The code is run for a given 
frequency and interblade phase angle for a given 
mode. The output is the real and imaginary 
components of unsteady pressures. This is repeated 
for NM modes of vibration. In the case of forced 
response for a gust, the code is run only once for the 
unsteady excitation. 

The mode shape values are interpolated onto the 
aerogrid before running LINLFUX as follows: At 
each CFD grid point on the blade surface, the 
distance to the nearest three finite-element grid points 
is calculated. Then, the modal deflections at these 
three nearest neighbors are used in a bi -linear 
interpolation scheme to calculate the interpolated 
value of the modal deflection at that CFD grid point. 
This is done considering the blade undeflected 
position as the reference position. The interpolated 

modal deflections <5 are stored and are used by 
LINFLUX. 

Calculation of elements of {AZ) r } 

The elements of {AD} are due to the external 


unsteady disturbances coming on to the blade. These 
may be due to entropic, vortical and wake 
interactions. These disturbances are independent of 
the blade vibratory motion, and therefore contribute 
only to the forced response aspect of the aeroelastic 
problem. The LINFLUX solver can calculate the 
unsteady pressures due to these excitations for a 
given interblade phase angle, r. Let p be the unsteady 
pressure due to these excitations on the blade surface. 
Then the elements of {AZ) r } are given by 

{AD r } - j P r dA.5 i (12) 

where 5- is the i* mode shape. 

Calculation of elements of [A.] 

The elements of M are due to the vibration of the 

blade in a given mode. It requires moving the 
aerogrid in the given mode in the aerodynamic 
calculations, and a given frequency, and calculating 
the unsteady aerodynamic forces. The LINFLUX 
solver can accept the modal values at each grid point 
and can modify the grid as per the modal values. 
Once the unsteady aerodynamic pressures are 
calculated, the elements of {A} are given as 

A rij =\5 i .dAp rj (13) 

where p r} is the unsteady pressure due to blade 
vibration in j lh mode for the rth interblade phase 
angle, S i is the i lh modal deflection, and dA is the 
elemental area. 

When running LINFLUX the pressure is usually 
output as a non-dimensional quantity, which is 
multiplied by (p*Cl ott ) to get the actual pressure 
values, where p is the reference air density and a ^ 
is the reference speed of sound. In calculating area 
and modal values, attention should be given to the 
non-dimensionlization to length scaling. If L is the 
reference length used for scaling the geometry and 
the modal values, the final Ay have to be multiplied 
by L**3 ( L for modal values and L**2 for area). It 
should also be noted that the modal values (mode 
shape) when input to LINFLUX should be divided by 
the reference length. 

Stability calculation 
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Aerodynamic work approach 

A work- per-cycle approach can be used to determine 
aeroelastic (flutter) stability. Using this approach, the 
motion of the blade is prescribed to be a harmonic 
vibration in a normal mode with a specified 
frequency. The vibration frequency is typically the 
natural frequency for the mode of interest, but some 
other frequency can also be used. The aerodynamic 
forces acting on the vibrating blade and the work 
done by these forces on the vibrating blade during a 
cycle of vibration are calculated. If work is being 
done an the blade by the aerodynamic forces, the 
blade is dynamically unstable, since it will result in 
extraction of energy from the flow, leading to an 
increase in amplitude of oscillation of the blade. 
Note that coupled mode flutter cannot be modeled 
with this approach. 

To determine aeroelastic stability using the work-per- 
cycle approach, the blade motion is specified to be 
harmonic: 

q(t) = q 0 sm(COt) (14) 

where q o is the amplitude of motion and CO is the 
vibration frequency. 

The work-per-cycle, W, done on the blade is 
calculated as: 

W = j>J pdA.(dX / dt)dt (15.1) 

s 

or, 

W = |J pdA.5q 0 co cos (QX)dt ( 1 5.2) 

S 


stable, but when W > 0, energy is gained by the blade 
and it is unstable. 

Eigenvalue approach 

An eigenvalue approach allows one to investigate 
coupled mode flutter. The flutter frequencies and 
flutter modes are calculated rather than assumed as in 
aerodynamic work approach. For a stability 
calculation (flutter), the motion-independent forces 
{■40,} are set to zero, and Eq. (10) can be written as 

■ ® j [m] {r} . M{r}- [ a,] m={o} a® 

Dividing Eq. 16 with an assumed frequency, C0 o 2 

-{co/co 0 f [Af]{y} + [ [AT]-[A r ]]{K}/ co 2 


= {0} (17) 

Rearranging, the equations can be written in the 
standard eigenvalue problem as: 

[[nhmpr} = { 0 } 00 

where 

[ P r] = [[K]-[A r ]]/(0 o 2 (19) 

[Or] = W (20) 

and 

y = (6)/fit) 0 ) 2 (21) 


The solution of the above eigenvalue problem results 
in NM complex eigenvalues of the form 


where, p = p(x,y y z,t) is the unsteady pressure on the 

blade surface due to blade vibration, A is the blade 
surface area vector pointing into the blade surface, 

J is the integral over the blade surface, ^ is the 
s 

integral over one cycle of blade vibration. The 
equation for work per cycle using the linearized 
unsteady aerodynamic equations is given in Ref. 18. 


. CO _ 

CO 

o 

The real part of the eigenvalue (fl ) represents the 
damping ratio, and the imaginary part ( V ) 
represents the damped frequency; flutter occurs if 

H — 0 for any of the eigenvalues. 




( 22 ) 


The work-per-cycle is an indicator of aeroelastic 
stability. The blade is dynamically unstable if the 
work done qh the blade during a cycle of blade 
vibration is positive. In other words, when W < 0 
energy is dissipated by the blade and the system is 


For the tuned cascade, the stability of each phase 
angle mode is examined separately. The interblade 
phase angle is fixed at one of the values given by 
equation (3), and the NM x NM eigenvalue problem 
is solved. The value of interblade phase angle is then 
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changed, and the procedure is repeated for each of the 
N permissible values. The critical phase angle is 
identified as the one, which results in the lowest 
flutter speed. 

Response Calculation : 

The aeroelastic response of the blades induced by 
wakes is calculated from equation (10) as 

W-[[P,]-mV{AD,} ( 23 ) 

The amplitude of the blade, (X), is obtained by 
summing contributions of { Y) from all the modes. 

Results 

Results are presented for selected three-dimensional 
geometry’s in this section. Results presented here are 
meant to demonstrate the state of development of the 
code. The procedure to do a complete aeroelastic 
analysis is shown in Fig.l. The procedure consists of 
running six codes, (I) a free vibration analysis and the 
corresponding processing code, RDVIB, (II) 
TURBO-AE, (III) INTERFACE, (IV) a PRE- 
processor, (V) LINFLUX, and (VI) a POST- 
processor. The function of each program is explained 
below in six steps and shown in table 1. In step 1 a 
vibration analysis is done for the blade. The output is 
structural gird coordinates, generalized masses, 
natural frequencies, and mode shapes. RDVIB reads 
this output and writes in a format required by the pre- 
processor. In step 2, a steady aerodynamic solution is 
obtained from TURBO-AE for the blade. The steady 
aerodynamic solution is written as a database. Step 3 
consists of running the INTERFACE program to 
rewrite the steady data base to the format required by 
LINFLUX. In step 4 the pre-processor is run. This 
will interpolate modal values at the structural grid 
onto the aerogrid, for all the modes. In step 5, 
LINFLUX is executed for the mode of interest and 
frequency. At this stage LINFLUX also calculates 
work per cycle for the mode of interest. In step 6 the 
post processor is run, to calculate the generalized 
forces, to calculate flutter using the eigenvalue 
approach and to calculate forced response 
amplitudes. 

Two examples are given below to show that the 
whole process is working as intended. In the present 
paper, the verification of unsteady pressures, and 
eigenvalue calculations is given. The unsteady 
pressures calculated using LINFLUX code are 


compared with those obtained from TURBO-AE and 
LINSUB codes. TURBO-AE is a three-dimensional 
Euler aeroelastic code developed in Ref.23, and 
LINSUB is a two-dimensional unsteady cascade 
aerodynamic code based on linear unsteady 
aerodynamic theory of Smith, Ref. 3. The eigenvalue 
calculations are compared with those obtained from 
ASTROP2 code of Ref. 8. The ASTROP2 code uses 
strip theory to integrate two-dimensional 
aerodynamic forces on to a three dimensional 
structure. The aerodynamic forces are calculated 
about leading edge, using linear unsteady 
aerodynamic theory of Ref. 4. It uses eigenvalue 
approach to calculate flutter stability. The 
verification of forced response amplitudes will be 
presented in a subsequent paper. 

Helical Fan 

The helical fan configuration consists of a rotor with 
twisted flat plate blades enclosed in a cylindrical duct 
with no tip gap. This configuration was developed 
by researchers to provide a relatively simple test case 
for comparison with two-dimensional analyses. Note 
that there is no experimental data available for this 
configuration. 

The parameters for this configuration are such that 
the mid-span location corresponds to a flat plate 
cascade with a stagger angle of 45 degrees, unit gap- 
to-chord ratio, operating in a uniform mean flow at a 
Mach number of 0.7 parallel to the blades. The rotor 
has 24 blades with a hub/tip ratio of 0.8. The radius 
at the hub is 8.619 cm (3.395 inches ) and the radius 
at the tip is 10.775 cm (4.244 inches). The inlet flow 
(axial) Mach number used in this calculation is 0.495, 
and the rotation speed of the fan is 16,962.4 rpm, 
which results in a relative Mach number of 
approximately 0.7 at the mid span section. 

The grid used for the calculations is 141x11x41 in 
each blade passage. On each blade surface, 81 points 
are located in the chordwise direction, 1 1 points in 
the span wise direction, and 41 points in blade to 
blade direction. The inlet and exit boundaries are 
located at an axial distance of approximately 0.7 
chord lengths from the blade leading and trailing 
edges. A steady solution is obtained first with 
TURBO-AE solver (step 2). Program INTERFACE 
is run to convert the steady aerodynamic solution 
data to the form required by LINFLUX (step 3). 
Next unsteady pressures are calculated for harmonic 
blade vibration in plunging and pitching modes. 
Mode shapes (step 1) are prescribed ( i.e. no 
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structural analysis is done ) such that the amplitude of 
vibration does not vary along the span. This choice 
of mode shapes is meant to reduce the three 
dimensionality of the unsteady flow field for ease of 
comparison with two-dimensional analyses. 

LINFLUX can be run by an input choice, iflut=l, by 
having these prescribed mode shapes internally 
calculated. However, the purpose of this effort is to 
run for arbitrary mode shapes that are calculated 
outside LINFLUX and to see that they are correctly 
read by LINFLUX. To check this, the mode shapes 
are obtained by running LINFLUX with iflut=l and 
are written on unit 3. PRE is then run to interpolate 
the mode shapes onto the aerogrid. But for this case 
this will not change the modal values at the aerogrid 
points since the modal values are obtained form 
LINFLUX to start with. So actually step 4 is not 
required for this example. After running PRE the 
interpolated mode shapes are written on unit 27. 
Now LINFLUX is run with iflut =2, which reads the 
externally supplied mode shape values at aero-grid 
points. LINFLUX is run for unsteady pressures and 
compared with those obtained with iflut=l and 
published results. 

Figures 2 and 3 show the unsteady pressures at mid 
span obtained for a reduced frequency, G), of 
1.0(0) = 0) o L/U = C0 o L/a oo M r = W/M r , 

where CO = (0 0 Li a w , and L is the reference length. 
It should be noted that CO is the input to LINFLUX.). 
For the present calculations the value of speed of 
sound, is assumed as 13707 in/sec, reference 
length is the chord which is 1.0 inch giving, a value 
of 1527 cycles/sec for the assumed frequency C0 o . 

Figure 2 is for plunging motion and Fig 3 is for 
pitching motion about leading edge. Fig 2.1 is for 
zero phase angle, fig. 2.2 is for 180 degrees phase 
angle. Figures for 3.1 and 3.2 show the unsteady 
pressure distributions for pitching motion for zero 
and 90 degrees phase angle respectively. The 
pressures are compared with those obtained from 
LINSUB code based on Smith theory. Ref. 3. To 
compare with LINSUB results, the unsteady pressure 
output from LINFLUX is divided by square of the 
Mach number times the amplitude of oscillation. The 
pressure distribution predictions from LINFLUX 
show excellent agreement with those from LINSUB. 

Next the unsteady pressure distributions from 
LINFLUX are compared with those obtained from 
TURBO-AE in figures 4 and 5. Figure 4 is for 
plunging motion for zero and 180 degree phase 


angles and figure 5 is for pitching motion for zero 
and 90 degrees phase angles. Again, excellent 
agreement is seen between the predictions. 

Table 2 shows the work-per-cycle calculations for 
this case, and comparison with TURBO-AE, Ref. 23. 
The work per cycle for plunging motion shows 
excellent agreement, however, the calculations for 
pitching motion show a lot of difference. This may 
be due to an error in normalization, since unsteady 
pressures shown earlier agreed well for this case. 
Some of these results were also reported in Ref. 18. 
They are included here for completeness, and to show 
that the present authors successfully ran LINFLUX 
code as intended, and for externally calculated mode 
shapes. 

Next, the calculated unsteady pressures are used by 
POST to calculate generalized forces and 
eigenvalues. The generalized forces are calculated 
using the prescribed mode shape used for the 
unsteady aerodynamic pressure calculation. To 
validate the eigenvalue calculations, the present 
results are compared with those obtained from the 
ASTROP2 code. As mentioned before, the 
ASTROP2 code uses strip theory to integrate two- 
dimensional aerodynamic forces on to a three 
dimensional structure. The aerodynamic forces are 
calculated about leading edge, using theory of Ref. 4. 
It uses eigenvalue approach to calculate stability. As 
unsteady pressures are already compared above with 
Smith theory, this aspect will give validation to the 
calculation of generalized forces and eigenvalues. 
Table 3 shows the eigenvalues obtained by the 
present post processor and the ASTROP2 code for an 
assumed frequency of CO =0.7 and for a = 0 
degrees. They show excellent correlation. As 
expected, the eigenvalues indicate that the blades are 
stable, which was also indicted by the work per cycle 
approach in TURBO-AE and LINFLUX programs. 
This validates the routines in the post-processor. 

Energy Efficient Engine , E3 

In this section, an example is considered where the 
mode shapes are not prescribed as for helical fan but 
obtained from a modal analysis. The configuration 
selected is derived from the Energy Efficient Engine 
(E-cubed) fan rotor. The E-cubed program was 
established by GE Aircraft Engines under NASA 
sponsorship in the 1980’s to demonstrate component 
technologies necessary to achieve higher efficiencies 
and reduce environmental effects in future subsonic 
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turbofan engines. Details regarding design and 
performance tests have been presented in Ref. 25, 

The fan rotor has 32 blades with a tip diameter, D, of 
210,8 cm (83 inches). The inlet Mach number is 0.5. 
The fan is rotating at a speed of 2727 rpm. The CFD 
grid is a simple sheared grid of size 43x15x1 1 with 
15 points on the blade surface in both the chordwise 
and spanwise directions. It is assumed that the tip 
gap is zero. The finite-element structural dynamics 
data is for a structural grid with 224 points. 

First a structural dynamic solution is obtained giving 
the modal frequencies and mode shapes, (step 1). The 
first mode is predominately bending with a frequency 
of 72.4 Hz, and the second mode is predominately 
torsion with a frequency of 163.8 Hz. 

A steady solution is obtained for this configuration 
using TURBO- AE. (step 2). The steady data base is 
processed through INTERFACE to get formatted 
data for LINFLUX (step 3). Then PRE is run to 
interpolate the modal values onto the aerogrid (step 

4) . The interpolated modal values are verified by 
plotting the deflected blade shape (not included here) 
and by manually checking the tip displacement 
values. Overall, the interpolated deflections matched 
the original deflections except for slight differences 
in the blade tip region. 

The unsteady aerodynamic solutions are then 
performed using LINFLUX. LINFLUX is run for a 
frequency of oscillation equal to the first mode (step 

5) . The non-dimensional frequency CO^CoD/a^ 
input to LINFLUX is 2.732. Figures 6.1 through 6.4 
show the unsteady pressures plotted for mid span for 
zero, 90, 1 80 and 270 degree phase angles. They are 
compared with TURBO-AE results. The TURBO- 
AE results are obtained using 100 steps per cycle. It 
can be seen that there is considerable difference 
between the two results, even though the trends are 
same. The difference is attributed to the fact that 
LINFLUX requires a very good steady solution, 
especially near leading edge, to get a comparable 
solution with TURBO-AE predictions. 

Figure 7.1 and 7.2 shows the unsteady pressures 
compared for the second mode. The non- 
dimensional frequency 65 = 6)D/a oo input to 
LINFLUX is 4.183. It is to be nofed that the value of 
the frequency parameter calculated with the second 
frequency is 6.18263. However, with this value the 
unsteady solution could not be obtained with 


LINFLUX. Therefore, the value is reduced to 4.183, 
where an unsteady solution could be obtained for a 
=0 and 270 degrees. Again, even for this value of 
CO , solution could not be obtained for a = 90° and 
1 80°. There may be two reasons for this behavior: (1) 
the steady solution is not sufficiently converged for 
LINFLUX run, (2) there is an inherent problem with 
LINFLUX that requires experience with CFL 
numbers, grid etc. 

Table 4 shows the work per cycle calculations from 
LINFLUX compared to TURBO-AE results. Both 
show a considerable difference. This is believed to 
be due to the coarse grid used in the calculations, and 
further calculations should be done with finer grid for 
this fan. Table 5 shows the eigenvalues calculated 
for this fan from POST and ASTROP2 for W = 
4.183 or C0 o = 111 Hz and for a = 0 degrees. There 

is some difference which is expected as the linear 
unsteady aerodynamic analysis in ASTROP2 neglects 
the effect of steady loading on flutter. 

Concl ud i ng Remarks 

The status of a turboamchinery aeroelastic code 
LINFLUX-AE is presented. The code is based on 
running four codes, one for steady aerodynamic 
solution, a linearized unsteady aerodynamic solution 
and for flutter and response calculations. Other 
modules required for connecting these solutions are 
developed and verified. These include a 
preprocessor to read structural vibration data, an 
interface for reading steady aerodynamic solution 
data and routine to interpolate mode shapes on to 
aero-grid, and routines for calculating eigenvalues 
and response amplitudes. 

The code is verified by running for two three- 
dimensional geometry fans. A helical fan with flat 
plate airfoils provided comparison with published 
results and a E-cubed fan providing calculation for a 
realistic fan. 

It was noted that for the helical fan with flat plate 
geometry, the unsteady pressures calculated from 
LINFLUX agreed well with LINSUB and TURBO- 
AE. However, the work per cycle calculations for 
pitching motion are not identical between LINFLUX 
and TURBO-AE. The flutter eigenvalues for an 
idealized flutter analysis showed the eigenvalues 
from LINFLUX-AE are close of those obtained from 
ASTROP2. 

For the E-cubed fan geometry, the unsteady pressures 
from LINFLUX and TURBO-AE showed 
considerable difference. This may be due to the fact 
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that the steady solution obtained from TURBO- AE is 
not a sufficiently accurate solution for LINFLUX 
code. This affected the work per cycle calculations 
also. In addition, unsteady solutions could not be 
obtained for some phase angles and some frequencies 
with LINFLUX. It is suggested that all the 
calculations for E-cubed fan be redone with a steady 
solution obtained with a fine grid. Further runs with 
some more 3D geometry’s and validation of forced 
response calculations are planned for future. 
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Table 1: LINFLUX-AE modules and their function 


steps 

program name 

function 

1 

RDVIB 

(1) do free vibration analysis. At present done outside this program 
through ANSYS, NASTRAN etc. 



(2) read vibration analysis output file and rotate the grid and mode 
shapes to aero-grid coordinate system, write on UNIT 3 

2 

TURBO-AE 

do steady aerodynamic analysis 

3 

INTERFACE 

convert TURBO-AE output to data format required by LINFLUX, 
write on UNIT 51 

4 

PRE 

interpolate (calculate) mode shape values at each aero-grid point on 
both the airfoil surfaces for all modes (UNIT 27); 
out.mode. 1 ; out.mode.2 out.mode.NM 

5 

LINFLUX 

calculate unsteady pressures for the given mode, frequency and 
interblade phase angle; repeat for all modes, frequencies and 
interblade phase angles (UNIT 94); 

out.prss.l and out.prsp.l; out.prss.2, out.prsp.2 

6 

POST 

calculate generalized forces; flutter eigenvalues and response 
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Table 2: Comparison of work/cycle, CO = 0.7, Helical fan, M=0.5, Q =16962 rpm 


Plunging 


sigma 

ND 

LINFLUX 

TURBO-AE 

0 

0 

-1.6854-02 

-1.6907-02 

90 

6 

-2.5382-02 

-2.5351-02 

180 

12 

-4.616-02 

-4.5597-02 

270 

-6 

-2.863-02 

-3.5906-02 


Pitching about mid-chord 


sigma 

ND 

LINFLUX 

TURBO-AE 

0 

0 

-0.13928-03 

-0.26938-03 

90 

6 

-0.13334-03 

-0.24221-03 

180 

12 

-0.33714-03 

-0.06438-02 

270 

-6 

-0.35759-03 

-0.07989-02 


Table 3: Eigenvalue comparison, Helical fan 


mode 

LINFLUX-AE 

ASTROP2 


Eigenvalue 

Eigenvalue 

Hi 

real 

imaginary 

real 

imaginary 

m 

-0.6533-06 

0.1527+04 

-0.4637-06 

0.152704+04 

2 

-0.1024-08 

0.4500+04 

-0.8846-08 

0.45000+04 


Table 4: Comparison of work/cycle, E-cubed fan, M=0.5, Q= 2727 rpm 


mode 1 ; CO = 2.732 


sigma 

ND 

LINFLUX 

TURBO-AE 

0 

0 

-0.18668-03 

-0.7225-03 

90 

8 

-0.87189-03 

-0.3526-02 

180 

16 

-0.2394-02 

-0.6938-02 

270 

-8 

-0.18337-02 

-0.3836-02 


mode 2; CO =4.183 


sigma 

ND 

LINFLUX 

TURBO-AE 

0 

0 

-0.7853-03 

-0.4225-02 

90 

8 

* 

0.2108-03 

180 

16 

* 

-0.5518-02 

270 : 

-8 

-0.53142-03 

-1.1175-02 


* the solution did not converge; 
got into numerically unstable region 


Table 5: Eigenvalue comparison, E-cubed fan 


mode 

LINFLUX-AE 

ASTROP2 


Eigenvalue 

Eigenvalue 


real 

imaginary 

real 

imaginary 

1 

-0.1001-00 

0.7217+02 

-0.5242+00 

0.7191+02 

2 

-0.7807-01 

0.1632+03 

-0.3021+01 

0.1635+03 
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Figure 1 : LINFLUX-AE flow chart 
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Figure 2.1 Unsteady pressures for plunging motion, 
a =0°, helical fan, G) =0.7, Mr=0.7 



Figure 3.1 Unsteady pressure distribution for pitching 
about leading edge, a=0°, helical fan, O)=0.7, 
Mr=0.7 



0 0.2 0.4 0.6 0.8 1 

distance along chord 


Figure 2.2 Unsteady pressures for plunging motion, 
a =180°, helical fan, 5) =0.7, Mr=0.7 



distance along chord 


Figure 3.2 Unsteady pressure distribution for pitching 
about leading edge, c =90°, helical fan, <jj=0.7, 
Mr=0.7 
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Figure 4.1 Unsteady pressure comparison, helical fan, 
LINFLUX vs TURBO-AE, plunging, a =0°. 



Figure 5.1 Unsteady pressures, pitching about mid- 
chord, helical fan, LINFLUX vs TURBO-AE, cr =0° 
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Figure 4.2 Unsteady pressure, helical fan, 
sigma=180, LINFLUX vs TURBO-AE, o =180° 
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Figure 5.2 Unsteady pressures, LINFLUX vs 
TURBO-AE, helical fan, a =90°, pitching about mid- 
chord 
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